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Abstract 



One of the fundamental assertions of statistical mechanics is that the time average of a 
physical observable is equivalent to the average over phase space, with microcanonical 
measure. A system for which this is true is said to be ergodic and dynamical properties 
can be calculated from static phase-space averages. Dynamics of a system which is fully 
integrable, that is has as many conserved quantities as degrees of freedom, is constrained 
to a reduced phase space and thus not ergodic, although it may relax to a modified 
equilibrium. 

In this thesis, we present a comprehensive study of chaos and thermalization of the 
one-dimensional Bose-Hubbard Model (BHM) within the classical field approximation. 
This model describes the dynamics of quantum degenerate gases in a lattice for sufficient 
occupation of every momentum mode and weak two-body scattering, and is of interest 
because of experimental advances of cooling and trapping alkali atoms in the quantum 
degenerate regime. 

We study chaos and its relation to thermalization. Two quantitative measures are 
compared: the ensemble-averaged Finite-time Maximal Lyapunov exponent, a mea- 
sures of chaos and the normalized spectral entropy, a measure of the distance between 
the numerical time-averaged momentum distribution and the one predicted by thermo- 
dynamics. A threshold for chaos is found, which depends on two parameters, the nonlin- 
earity and the total energy-per-particle. Below the threshold, the dynamics are regular, 

x 



while far above the threshold, complete thermalization is observed, as measured by the 
normalized spectral entropy. 

We study individual resonances in the Bose-Hubbard model to determine the cri- 
terion for chaos. The criterion based on Chirikov's method of overlapping resonances 
diverges in the thermodynamic limit, in contrast to the criterion parameters inferred from 
numerical calculations, signifying the failure of the standard Chirikov's approach. 

The Ablowitz-Ladik lattice is one of several integrable models that are close to the 
BHM. We outline the method of Inverse Scattering Transform and generate the integrals 
of motion of the Ablowitz-Ladik lattice. Furthermore, we discuss the possible role of 
these quantities in the relaxation dynamics of the BHM. 
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Chapter 1 



Introduction 

1.1 Introduction to Dynamical Systems 

The study of thermalization in nonlinear systems dates back to the early numerical stud- 
ies of coupled anharmonic oscillators by Fermi, Pasta, and Ulam (FPU). At the time 
it was expected that any amount of nonlinearity, no matter how small, in a system 
with many degrees of freedom would lead to thermal behavior. Quite unexpectedly, 
thermalization was not observed. Later work found a threshold in the nonlinear cou- 
pling strength, with energy equipartition occurring above the threshold. The absence 
of thermalization for small nonlinearities has been explained in terms of the presence 
of a threshold for the onset of widespread chaos determined by Chirikov's criterion of 
overlapping resonances and also in terms of the closeness to a fully integrable model, 
the Korteweg-de-Vries equation. Additional studies on thermalization and approach to 
equilibrium have been carried out in many classical field theories. 

Fermi, Pasta, Ulam (FPU) and Tsingou performed one of the first numerical exper- 
iments in the 1950's (Fermi et al, 1974). Their discovery led to significant work and 
important discoveries in the field of nonlinear dynamics. The motivation for the experi- 
ment was to study the approach to equilibrium of a nonlinear system. The contemporary 
expectation was that an system with a large number of degrees of freedom would exhibit 
thermal behavior in the presence of any non-linearity, no matter how small. In particular 
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they were looking to observe the Fourier heat law for conduction. They numerical inte- 
grated a system of one-dimensional anharmonically coupled oscillators with quadratic 
and cubic forces. 

The results found by FPU were completely suprising. In order to understand the 
prevalent view at the time and the significance of their results, let us take some time to 
review important concepts in Hamiltonian systems and statistical mechanics. 



1.2 Dynamical Systems 

Given a system whose state is completely known at some time to, what can be said 
about the state at some future time t > t(p. What about its state at some time in the past, 
t < *o? Classical mechanics addresses these questions. We begin with a short review of 
Hamiltonian dynamics, referring to Tabor (1989). 

Hamiltonian Systems Consider a generic Hamiltonian of an N-particle system, 
H{{quPi}) where are the coordinates of the N-particles and {pt} are the corre- 
sponding momentum. The equations of motion are given by 

opt dqi 

These equations, called the Hamilton's equations of motion, uniquely determine the 
time evolution of the state, from the initial state, which is specified by a complete set 
of values {quPi}. Note that canonical coordinates and momentum satisfy Liouville's 
theorem, so that a volume element in phase space moves like an incompressible fluid 
under the Hamiltonian flow. 
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Poisson Brackets For any dynamical quantities f,g in a system with canonical coor- 
dinates {quPi}, the Poisson brackets are defined as 



if i-y^^ d/ jjg 

~~ y dpi dqi dq t dpi ' 



From this one can write the time dependence of a function /(/?, g, t) as 



dt i \dqi dt dpi dt J dt 
y\dqidpi dpidqt) dt 



Thus any quantity that does not explicity depend on time is a constant of motion if the 
Poisson bracket of that quantity vanishes. The Poisson brackets are antisymmetric and 
satisfy the identity 

{f,gh} = g{fM + {Lg}h- 

Canonical Transformations Oftentimes it is more convenient to work in one coor- 
dinate system than another. A canonical transformation is one in which the canonical 
form of Hamilton's equations is preserved. The phase volume is preserved under a 
canonical transformation. Thus the Jacobian, must be unity. Consider a transforma- 
tion from one set of phase-space variables {puqi) to a new set of variables (PuQi) the 
preservation of phase volume is expressed as 

J dp\dqi j dp 2 dq 2 ... j dp N dq N = J dP\dQ\ j dP 2 dQ 2 ... j dP N dQ N 
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and the Jacobian must satisfy the condition 

d(pi,...,p N ,qi,...,q N ) d(Pi,...,P N ,Qi,...,Q N ) 



B(Pi,...,Pn,Qi,—,Qn) d(pi,...,PN,qi,--,qN) 



= 1 



(Tabor, 1989). 

One of the representations of the canonical transformation, is through the type 2 gen- 
erating function 4>(0,/), which transforms 7,0 — >■ /,0. The Hamiltonian is transformed 
according to 

H(I,Q) =H(I,Q) + d -^ 
3<l> 

1 = 30 (L2) 

a* 
9 ~W 

1.2.1 Integrable Systems 

An integrable system has as many independent constants of the motion as degrees of 
freedom. Stated another way, a systems with n degrees of freedom is completely inte- 
grable is there exist n integrals of motion which are in involution. 

{Fi, H} = 0, i = 1 • • • n, {Fj} are constants of the motion 

{Fi, Fj} = 0, i = 1 • • • n, j = 1 • • • n {Fi} are in involution 

Ford makes the distinction that these must be in fact well-behaved integrals of 
motion (Ford, 1992). All one dimensional systems are integrable (Tabor, 1989). Exam- 
ples of integrable systems with many degrees of freedom include the Toda lattice 
(infinite-dimensional countable), the Korteweg-de-Vries equation, sine-Gordon model 
and continuous nonlinear Schrodinger equation (infinite-dimensional continuous). 
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1.2.2 Chaos 

The term chaos is commonly used to describe systems where the motion appears ran- 
dom, erratic or unpredictable. It is sometimes called "deterministic chaos" to emphasize 
that it refers to dynamical systems in which the motion is governed by deterministic 
equations of motion. A key feature of chaotic systems is extreme sensitivity to initial 
conditions. Imagine two identical systems that are governed by the same equations of 
motion and with initial conditions that are only slightly different. In a chaotic system, 
the initial different will grow rapidly and after some time the states of the two systems 
will be entirely different from one another. In contrast, in a regular system, the differ- 
ence will grow slowly, so that two remain highly correlated. The difference between 
chaotic and regular motion is not only quantitative, but qualitative. The very functional 
form of the divergence is different: chaotic trajectories diverge exponentially while reg- 
ular trajectories diverge linearly. 

Lyapunov Exponents A chaotic system is characterized by local instability, specifi- 
cally by exponential divergence of trajectories neighboring in phase space. This rate of 
divergence of is measured by the Lyapunov exponent. The maximal Lyapunov Exponent 
(MLE) of a system is defined as 

Mx„)=lim lim lln'^-'W" 
r^-xo^xo t ||x -Xo|| 

where x(t) and x(t) are two phase space trajectories (Tabor, 1989). Chaotic motion 
is characterized by positive MLE, X > 0, in which case trajectories that are initially 
close in phase space will diverge exponentially. On the other hand, a zero MLE, X = 0, 
indicates regular motion and linear divergence. The Lyapunov exponent is a function 
of an initial state and is thus a measure of local chaos. A system may have a mixed 
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phase space, consisting of both chaotic and regular regions (Zaslavsky, 1999). Thus 
initial states in the regular regions will have zero Lyapunov exponent while initial states 
in chaotic regions have a positive Lyapunov exponent. As a system transitions from 
regular to globally chaotic behavior, the volume of phase space corresponding to chaotic 
regions grows while the regular regions shrink. When the chaotic regions dominate the 
phase space, there can still exist regular regions, called islands of stability, which are 
surrounded by the chaotic sea. 

1.3 Statistical Mechanics 

Thermodynamics is a phenomenological theory that has very successfully described the 
equilibrium properties of isolated systems with many degrees of freedom. In this sec- 
tion we give an overview of basic concepts of statistical mechanics. The field arose in 
the study of macroscopic systems of atoms with N ~ 10 23 particles and correspondingly 
large volumes. One of the great success of thermodynamics is the Maxwell-Boltzmann 
distribution for ideal gases. It can be derived in two ways: (1) through the Boltzmann 
transport equation and (2) through the most probable distribution. Systems are described 
not by position and momentum of ever particle in the system, but by probability distri- 
bution functions from which macroscopic properties such as pressure and temperature 
and volume can be determined. The aim of statistical mechanics is to derive the laws of 
thermodynamics from molecular dynamics. 

For a system with N particles with canonical coordinates , #2, • • • , q^N and conjugate 
momenta qx.qi, ■■■■ ) q3N, the 6N dimensional space spanned by these coordinates is the 
phase space or T space. A representative point is a point in this T space that specifies 
the position and momentum of each of the N particles. The representative point will 
also be called a microstate. For each microstate, there is a corresponding macrostate 
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that specifies properties of the system such as density, temperature, and pressure. Many 
microstates and in fact, an infinite number, correspond to the same macrostate. 

The dynamics of the individual N particles is governed by the Hamiltonian, H{p, q) 
where (p,q) = (p\,pi, ... 7 p3N,qi,q2, ■■■,13n)- The equations of motion are given by the 
Hamilton's equations of motion, 

dH 

dpi 
_BH 

These equations become quickly intractable for systems with large numbers of degrees 
of freedom. Instead systems will be studies by various ensembles which give the distri- 
bution of representative points in phase space. 

1.3.1 Statistical Ensembles 

A statistical ensemble is a collection of microstates which correspond to the same 
macrostate. Geometrically the ensemble can be described by the distribution of rep- 
resentative points in the T space with the density distribution function p(p,q,t) which 
is defined so that 

p(p,q,t)d 3N pd 3N q 

is the number of representative points in phase space volume d 3N p d 3N q at time t. 
Liouville Theorem 




The interpretation of the Liouville theorem is that the distribution of representative 
points in phase space move like an incompressible fluid. 




7 



Ensemble averages The ensemble average of an observable is give by 



(0) 



fd 3N pd 3N qO(p,q)p(p,q,t) 
fd™pd™q P (p,q,t) 



The time dependence of comes from the time-dependence of p. 

Postulate of Equal a Priori probability For a system in thermodynamic equilibrium, 
it is equally likely to be in any microstate (of the same phase volume) that satisfies the 
macroscopic conditions of the system. 

For an isolated system, the distribution is described by the microcanonical ensem- 
ble, which corresponds to all microstates with the same energy, volume and number of 
particles. The density of representative points in phase space is given by 



The question arises, what ensemble describes a system that is not in isolation, but 
instead is in equilibrium with another, larger system? A system in contact with a heat 
reservoir such that the temperature, volume and number of particles are constant is 
described by the canonical ensemble. Is is also possible to have a system where par- 
ticles can be exchanged with a larger system. Such a system, that is in contact with a 
heat reservoir and a particle reservoir is described by the grand canonical ensemble, in 
which temperature, volume and chemical potential are kept constant. 




E <H(p,q) <E+A 



otherwise 
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The most probable value of an observable is given by the value of O that the most 
members of the ensemble have. The most probable value and ensemble average are 
close if the mean square fluctuation 



Ergodic theorem Under certain conditions, a representative point in phase space will 
pass arbitrarily close to any other point in the accessible phase space, if one waits a 
sufficiently long time. Following from the postulate of equal a priori probability, the 
time average of some physical observable is equal to the ensemble average over phase- 
space, that is 



Statistical mechanics does not specify whether a system is ergodic or not and there 
is no generic test for ergodicity. For an ergodic system, a single trajectory will uni- 
formly cover the phase space. An integrable system will not be ergodic in the full phase 
space due to the additional conserved quantities that act as constraints. Ergodicity says 
nothing about the time-scale involved in covering the entire phase space. For typical 
thermodynamic systems, the size of the phase space is immense. A stronger condition 
than ergodicity is mixing (Tabor, 1989). In the long time limit, the values of a macro- 
scopic observable in a system that is mixing will equal the ensemble average, without a 
need to time average over the trajectory as in an ergodic system. 



(o 2 ) - (Q) 2 

(O) 2 



< 1 



is small. 
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1.4 FPU Model and Results 



Let us now return to the numerical experiments of FPU. Recall that the expectation was 
that a nonlinear system with many degrees of freedom would exhibit thermal behavior 
for any non-linearity, no matter how small. The stystem studies was a chain of one- 
dimensional anharmonically coupled oscillators with quadratic and cubic forces. The 
Hamiltonians of these two models can be written as a sum of 

H = H + H h (1.3) 

where the integrable Hamiltonian 

N-l N-l 
n= 1 n= 1 

is weakly perturbed by the non-integrable Hamiltonian Hi, which for the a— model is 

#1 = T £ (Xn+l-Xn? (1-5) 
5 n=i 



and for the (3— model is 



#1 = -7 £ (x n +l -Xn) 4 '• (1-6) 
4 n=\ 



The displacement of particle n from equilibrium is x n and a and (3 are nonlinear interac- 
tion parameters. The resulting equations of motions for the a- model are 

X n = (x n +l -2x n +X n -i) +CC[(x M+ l - X n ) 2 - (x n - X n -{ ) 2 ] (1.7) 

and for the [3-model 

X n = (X n +1 -2x n +X n -i) +$[(x n+ \ - X n f - (x n - JC„_i) 3 ]. (1.8) 
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The problem can be analyzed in terms of normal modes, Ak(t), which are the Fourier 
modes of the displacement, x n {t), 



In terms of the normal modes, Hq is a sum of harmonic oscillators and H\ is a perturbing 
term that couples the oscillators. 



For the linear system, that is when a = (3 = 0, there is no interaction between nor- 
mal modes, so that modes that are initially populated remain populated and modes that 
are initially unpopulated will remain unpopulated. What happens when nonlinearity is 
added? Would the coupling between the modes cause the energy to spread from a single 
mode to all of the other modes in the system? For the systems studied, thermalization 
would be marked by equipartition of energy among all of the modes, A k E^ = A' k E' k for all 
k, k'. The expectation was that thermal behavior would be observed. In Fig. 1.1 the time 
dependence of the normal modes is plotted for N = 32 oscillators with cubic forces and 
(3 = 8. The initial condition is given by a sine wave and the velocity is zero. The ends 
are fixed. The model preserves symmetry so that the effective number of particles is 16 
and even modes have zero energy. As can be seen from the plot only a few modes are 
active in the dynamics and there are recurrences. Additionally the period of recurrence 
was found to decrease with increasing nonlinearity. Later work found a super period, 
where almost all of energy (99% ) returns to the initial modes after about 80 000 T\ . 
(Tuck and Menzel, 1972). 




(1.9) 



Ho 



k 
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f U THOUSANDS OF CYCLES 



Figure 1.1: Time dependence of normal modes for cubic forces with (3 = 8. (Fermi 
etal, 1974) 

The failure of FPU to thermalize was very puzzling when first discovered and it was 
some time before explanations came out to explain the observed phenomena. The first 
explanation has to do with closeness to integrable systems and the second has to do with 
a stochasticity threshold with respect to the linear system. 

1.4.1 KAM Theorem 

A theorem outlined by Komologorov and subsequently proved by Arnold and Moser 
provided one resolution to the apparent paradox of FPU (Kolmogorov, 1954; Moser, 
1962; Arnold, 1963). Consider an integrable Hamiltonian that is weakly perturbed: 
H = Ho(l) + eH\ (1, 0) where Hq is integrable with constants of motion I and frequencies 
(Hi = ^and Hi is periodic in the original angle variables, H\ (1, + 2tc) = H\ (1, 0). The 
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motion of the unperturbed Hamiltonian corresponds to motion on an n— dimensional 
torus. It is assumed that the Hamiltonian is analytic on the complex domain and that the 
unperturbed Hamiltonian is non-degenerate, 



det 



d 2 H 



dlidlj 



Consider a frequency vector of the unperturbed Hamiltonian, (0* that is incommen- 
surate ((0* • k 7^ for all integer ki). The corresponding motion of the unperturbed system 
is on the torus 7o((Q*). 

One statement of the KAM theorem is 

Theorem 1.4.1 KAM Theorem If Hi is small enough, then for almost all &*, there exists 
an invariant torus T(&*) of the perturbed system such that T(&*) is close to Tq(&*). 
(Arnold and Avez, 1968). 

Stated informally, KAM showed that if the perturbation is sufficiently small, then 
almost all of the tori of the unperturbed motion are preserved and the resulting motion 
is quasi-periodic. 

1.4.2 Solitons and the Korteweg-de-Vries Equation 

Soon after the proof of the KAM theorem, Zabusky and Kruskal discovered solitary 
wave solutions to the Korteweg-de-Vries (KdV) equation, which they termed "soli- 
tons"(Zabusky and Kuskal, 1965). Solitons are solitary wave that preserve their shape 
both under free propagation and after collisions. Zabusky and Kruskal also show that 
that the continuum limit of the |3— FPU model is close to the KdV equation. In KdV the 
speed of the solitons depends only on their amplitude. The KdV equation is given by 

u t + uu x + d 2 u xxx = (1-10) 
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and was first found as a description of the motion of shallow water waves. Starting with 




Figure 1.2: Space-time diagram of soliton trajectories, uq = cos(7Uc). 8 = 0.022. T R is 
the recurrence time (Zabusky and Kuskal, 1965). 

a cosine pulse, the negative slope regions of u steepen, and then oscillations develop 
on the steep front and grow in amplitude, and finally each solitary wave or "soliton" 
moves at a constant speed, which is proportional to the amplitude. The trajectories of 
interacting solitons are shown in Fig. 1.2. The solid (dashed) lines represent the odd- 
(even-)numbered solitons. From the trajectories, it is clear that when solitons interact, 
they emerge with the same speed and thus shape. The dotted lines represent interactions, 
during which the joint amplitude is less than the sum of the individual amplitudes due 
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to the non-linear interaction. At the recurrence time Tr, all of the solitons arrive at one 
point in space and almost reconstruct the initial state. 

KdV was later shown to be fully integrable by the method of Inverse Scattering 
Transform (Gardner et al, 1967). The discovery of the solitons of KdV provides one 
route for explaining the absence of thermal behavior in FPU. The KAM theorem pro- 
vides another explanation. The nonlinearities of the FPU studies were insufficient to 
break the KAM tori and thus the motion remained quasi-periodic. Indeed later work 
showed that for larger nonlinearities, FPU did exhibit thermal behavior 

1.4.3 Chirikov Criterion 

A method for predicting the criterion for the onset of chaos in the FPU experiments was 
the theory of overlapping resonances developed by Chirikov in the context of plasma 
dynamics, and later applied to FPU (Izrailev and Chirikov, 1966). 

Consider a one-dimensional non-linear oscillator perturbed by an external periodic 
force. Write the unperturbed system in action-angle variables (7, 0) and the external 
field as a Fourier series. 

H = H (I)+£^V mn (I)e l ^ e+ '^ 

where Hq is the unperturbed Hamiltonian with frequencies Cfl(7) = -^f-, (|) = Clt + (|)o is 
the phase of the external force and = cof + (f)'. Resonances occur for the set of I r such 
that a>(I r )l = CLk. The main ingredients are then to: 

1. Assume that when the system is near a resonance, the resonance dominates the 
motion. Thus resonance are studied in isolation. 
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2. Make a canonical transformation from (7,0) — > (J,i\f) into the rotating reference 
frame of the resonance. (\\f measures deviations from resonance). The old and 
new coordinates are related by 

7=l(/-/ r ), y = lQ-k§. 

3. Integrate out the fast phase § motion. 

4. Expand the Hamiltonian about I = I r and keep terms up to second order in Hq(I) 
and zeroth order in V nm (I). 

The resulting Hamiltonian of this process is, 

J 2 

Hr= 2M +eV, '-* C0S M 

which is the Hamiltonian of a simple pendulum, with "mass" M = I 2 (^Jr j • In 
Fig. 1.3 a schematic is plotted for 7, \\f and \|/ coordinates generated by two different 
canonical transformations corresponding to two different resonances. The separatricies, 
which separate bounded and unbounded motion, are given by the black contour lines. 
An initial state inside the separatrix is captured by the resonance and both the action vari- 
able and the angle are bounded in phase-space. Outside of the separatrix, the angle is 
unbounded. When the resonances are well-separated in phase-space, that is f — J 3> A/, 
the motion near an individual resonance is dominated by that resonance. As the strength 
of the external drive increases, the width of the separatrix, A7, grows. Eventually the 
width of the two separatricies becomes comparable to the distance between the reso- 
nances. 
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Figure 1.3: Schematic of Chirikov's criterion of the onset of chaos for the kicked rotor. 

and J' ,\|/ are coordinates generated by two different canonical transformations. 
AJ and A/' are width of the separatricies. J — J' represents the distance between the 
resonances 

For two neighboring resonances, with resonant values Ia and /g, the "distance" 
between the resonances is given by Ib — Ia- For separatrix width AI, the ratio of these 
two quantities, 

separatrix width AI 



distance between resonances I B — Ia 
gives the condition for onset of chaos in the system 

Crossing the threshold corresponds to the overlap of the separatrices of the two reso- 
nances in phase space, so that the action variable is free to travel between resonances 
and thus explore all of phase space. Complete chaos emerges when all of the neighbor- 
ing resonances become coupled. 
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Chirikov applied this criterion to the FPU system. For a chain of N oscillators, and 
an excitation of momentum mode k, the conditions for chaos are given by: 



Low modes 



k<^N: 




(1.11) 



High modes (k « AO 



N-k<^N: 




(1.12) 



This criterion predicts that for initial excitation of low modes, there is only have stochas- 
ticity for large perturbations, while for the high mode have stochasticity even for small 
non-linearities when N is large. Data from FPU experiments show a threshold for ther- 
malization that is close to the prediction of Chirikov. 



1.5 Relationship between Chaos and Thermalization 



Statistical mechanics dictates the of equilibrium state of a system, but it does not tell 
us whether a system will thermalize or not. How is it that statistical behavior can arise 
from dynamics? 

Consider a system in equilibrium, with some constraint. Lift the constraint and let it 
evolve. What will the equilibrium state be, if in fact it reaches equilibrium? According 
to the second law of thermodynamics, when a constraint is lifted, the system moves 
to a state of greater entropy. To explain this further, consider a simple example a box 
of volume V with N particles, initially confined to 1/2 of the box. When Boltzmann 
first derived statistical mechanics for an ideal gas through kinetic equations, there were 
several objections that were raised (Zaslavsky, 1999): 
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Objection 1: Zermelo's Paradox (Recurrence) Poincare's recurrence theorem states 
that after sufficiently long times, any trajectory will pass arbitrarily close to any phase 
space point, including the initial state. This would contradict the expected increase in 
entropy. 

Objection 2: Loschmidt's Paradox (Reversibility) The equations of motion are 
invariant under time-reserval. If one simply reverses the velocities of particles, one goes 
back to the original state, a process that decreases entropy. 

Microscopic Origins of Macroscopic Irreversibility The second objection to the 
can be restated, how is it that the microscopic equations of motion are reversible, but the 
macroscopic behavior is irreversible? Lebowitz argues that this question was satisfac- 
torily settled years ago by Thomson, Maxwell and Boltzmann (Lebowitz, 1999). The 
essential ingredients of understanding this, according to Lebowitz are 

1. the vast difference in scales between microstates and macrostates 

2. initial conditions are special 

3. significance of probabilities 

Consider an isolated classical system of N particles. Let X be the microstate that 
completely specifies the system, X = (ri , pi , r2,P2, • • • , r#, Pat) in phase space T. Let M 
represent the macrostate of a system and let Ym be the the region of the phase space T 
that correspond to macrostate M. Note that there are many microstates X that correspond 
to the same macrostate M. The number of microstates that correspond to is so large as 
to make certain macrostates extremely unlikely. Consider the example of the gas of 
particles initially confined to half of a box. The initial state is indeed a special state. 
Once the constraint is lifted, the volume in phase space corresponding to all of the 
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particles in one-half of the box is vastly smaller than the volume of phase space that 
corresponds to a roughly equal distribution of the particles between the two halves of 
the box. Thus, the probability that it will return to the initial state is effectively zero. 
Boltzmann made a rough estimate of the Poincare recurrence times and found them to be 
much larger than the lifetime of the universe and thus irrelevant. This interpretation of 
the second law of thermodynamics explains why FPU expected to see thermal behavior 
in their system. However, it fails to account for the thermalization threshold found in 
FPU. Additionally, thermal behavior has been observed in systems with at few degrees 
of freedom, where this reasoning does not apply. 

Zavslasky argues that chaotic dynamics introduce mixing properties in a system 
(Zaslavsky, 1999) that well resolves these paradoxes. Furthermore, he calculates the 
distribution of Poincare recurrence times and concludes that they are irrelevant. While 
there is not a current consensus on the origins of statistical laws, this discussion high- 
lights some relevant questions, namely, 

• Is chaos necessary for thermalization? 

• What is the role of chaos in thermalization in systems with many degrees of free- 
dom? 

1.5.1 Thermalization in Classical Field Models 

Since the FPU studies on anharmonic oscillators, further studies on thermalization 
and approach to equilibrium have been carried out in several classical field theories, 
including recent studies on the classical (|) 4 model (Boyanovsky et al, 2004), nonlinear 
Klein-Gordon equation (NLKG) (Gerhardt et al. , 2002), nonlinear Schrodinger equation 
(NLSE) (Villain and Lewenstein, 2000; Herbst and Ablowitz, 1989), discrete nonlin- 
ear Schrodinger equation (DNLS) (Herbst and Ablowitz, 1989; Ablowitz et al, 1993) 
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equivalent to the Bose-Hubbard model, and Integrable Discrete Non-Linear Schrodinger 
equation (IDNLS)(Herbst and Ablowitz, 1989). 

No conventional thermalization is expected in the NLSE and IDNLS, which are both 
integrable. In NLKG, like in FPU, the ability of the system to reach thermal equilibrium 
in the course of time evolution emerges only when the degree of nonlinearity exceeds a 
certain critical value (see (Izrailev and Chirikov, 1966; Livi et al, 1985) for the thermal- 
ization threshold in FPU). On the contrary, the (|) 4 model eventually reaches equilibrium 
regardless of how small the nonlinearity is. 

There are several other studies on thermalization and chaos in system with a large 
number of degrees of freedom that are highly relevant to the work presented here. Livi 
et al. (1985) investigated the equipartition threshold in the FPU |3 model in the thermo- 
dynamic limit. For N oscillators, 64 < N < 512, the thermodynamic limit is simulated 
by initially exciting a block of modes, An, such that ^ remains constant. The threshold 
for equipartition of energy is found to be independent of the number of degrees of free- 
dom with respect to the relevant control parameter, the energy density, with a critical 
value of £ c ~ 0.35. They also calculate the Asymptotic Reynolds number, R, given by 



which is a measure of ratio of strength of nonlinear to linear terms. Again, there is 
universal behavior with R c ~ 0.03, consistent with findings for energy density. There is 
evidence that the threshold energy is independent of the mode excited, when a narrow 
range of energies is initially excited. Note that very long equipartition times (t — > °°) are 
not ruled out and they conclude that the results are relevant for long, but finite times. It is 
significant to note that the result that the threshold for FPU remains in the thermodynam- 
ics limit contradict the predictions of Chirikov's criterion of overlapping resonances. 
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Another study by the same group focuses on the relationship between chaotic 
dynamics and statistical mechanics in two nonlinear Hamiltonian systems, the FPU 
model of nonlinearly coupled oscillators and coupled rotators (Livi et al, 1987). For 
both systems thermodynamic quantities are computed analytically using ensemble the- 
ory and compared with dynamical results from numerical simulations. For the FPU 
model, there is qualitative agreement between ensemble-averages and time-averages, 
independent of the stochasticity. That is the system is ergodic in both the chaotic and 
regular region. For the rotator model, there is good agreement between the ensemble- 
averages and time-averages at low temperatures but not at high temperatures where the 
system is strongly chaotic. This result is explained in terms of localized chaos. In con- 
clusion it is possible that a system (a) is not chaotic, but is ergodic for some physically 
"relevant" quantities and also (b) is chaotic, but some observable are not ergodic. This 
study highlights the open questions in the relationship between stochasticity and ther- 
malization, particularly in systems with many degrees of freedom. 



1.6 Quantum Degenerate Gases - Ultracold Atoms 

Advances in the cooling and trapping of alkali atoms into the quantum degenerate 
regime has led to an explosion of experimental and theoretical studies of ultracold atoms. 
In recent years Nobel prizes have been awarded for advances in laser cooling techniques 
(Phillips and Metcalf, 1982; Chu et al, 1985; Aspect et al, 1988) and the subsequent 
observation of Bose-Einstein condensation (BEC) (Davis et al, 1995; Anderson et al, 
1995). The manipulation of atoms by electric and magentic fields offers unprecedented 
control over parameters in the system and the ability to address fundamental questions 
in physics. Numerous proposals have been put forth for quantum simulators and appli- 
cations have arisen in precision measurements. Studies in ultracold atoms have led to 
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fruitful collaborations across fields, such as condensed matter and quantum informa- 
tion. One application of these advances is using matter-wave interferometers based on 
ultracold atomic systems for high precision sensing of accelerations and gravitational 
fields (Gustavson et al, 2000; Durfee et al, 2006; Weiss et al, 1994; Wicht et al, 
2002). Fundamental questions in physics related to out-of-equilibrium dynamics and 
thermalization in classical and quantum integrable systems have also been studied in 
one-dimensional ultracold atoms. Quasi-one-dimensional systems have been realized in 
optical lattices (Paredes et al, 2004; Kinoshita et al, 2006) and on atom chips, where 
BEC's have been created and manipulated (Esteve et al, 2006; Schumm et al, 2005; 
Wang et al, 2005). 

1.6.1 Bose-Einstein Condensation 

Advances in laser cooling and trapping led to the realization of Bose-Einstein conden- 
sation (BEC) in alkali atoms (Davis et al, 1995; Anderson et al, 1995; Bradley et al, 
1995). BEC is a phase of matter, first proposed by Bose (1924) and Einstein (1925), 
in which there is macroscopic occupation of a single quantum state. Seventy year later 
BEC was created in Rb 87 gas (Anderson et al, 1995), sodium (Davis et al, 1995) and 
Li 7 (Bradley et al, 1995). Since the initial experiments, BEC has been observed in 
twelve species of alkali atoms as well as in Bose molecules (Yukalov, 2009). BEC was 
created by confining and cooling atoms to microkelvin temperatures with a magneto- 
optical trap (MOT), followed by evaporative cooling to nanokelvin temperatures. In 
Fig. 1.4 the velocity distribution of rubidium atoms is show prior to and after conden- 
sation. The velocity distribution of rubidium atoms is measured by turning off the con- 
fining trap, allowing the atoms to expand and performing a time-of-flight measurement. 
The leftmost plot shows the velocity distribution just before condensation. The center 
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Figure 1.4: Velocity distribution of rubidium atoms. Left: prior to condensation. Center: 
just after condensation. Right: after further cooling (Anderson et al, 1995) 

plot is the velocity distribution just after condensation, where the sharp peak in veloc- 
ity distribution clearly indicates the presence of the condensate. In rightmost plot, the 
system has been cooled further such that most of the atoms are in the condensate. The 
presence of the condensate was confirmed by the anisotropic velocity distribution due 
to the magnetic trap in contrast with the isotropic, thermal velocity distribution. BEC's 
demonstrate long-range phase coherence, confirmed experimentally by the observation 
of interference between two independent condensates (Andrews et ah, 1997). 

Historically BEC was defined for a uniform, ideal gas as the macroscopic occupation 
of a single quantum state in the thermodynamic limit, 

lim ^>0 

A'/V^const 

where N is the total number of particles, V is the volume and No is the occupation num- 
ber of a single quantum state (Yukalov, 2009). The question arises as to to define BEC 
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in a non-uniform system such as in presence of trap. Penrose and Onsager propose that a 
condensate is present when the largest eigenvalue of the single-particle density matrix is 
extensive (Penrose and Onsager, 1956). The Penrose-Onsager scheme is more general 
than idea of off-diagonal long-range order and applicable to both uniform systems and 
trapped systems. There is no true condensate in finite or ID systems, which will be the 
focus of this work, however there can be a quasi-condensate in ID, when the coherence 
length is much larger than the de Broglie wavelength (Castin, 2004). 

1.6.2 Bosons in Optical Lattices 

The versatility offered by optical lattices allows one to control parameters such as the 
interaction strength, lattice spacing and the dimensionality of the system. In particular, 
one-dimensional systems in cold atoms have been realized in optical lattices by tight 
confinement in two dimensions. 

The dynamics of ultracold bosons in optical lattices can be described the Bose- 
Hubbard model (BHM). The Bose-Hubbard Hamiltonian (Jaksch et ah, 1998) is 

H = - J E + \uY,ni{hi- 1). 

To derive this Hamiltonian, one begins with the Hamiltonian for bosonic atoms in an 
external potential 

H = J dW(*) (-|^V 2 + Vo(x) + Vr(x))\|/(x) 

+ 14naji_ r j3 w t (x)¥ t (x)¥(x)v(x) 
l m J 

where V|/(x) is a boson field operator, Vq(x) is the potential of the optical lattice, VV(x) 
is the trapping potential and a s is the s-wave scattering length. For single atoms, energy 
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eigenfunctions are Bloch wave functions. If the energies involved are much less than 
the excitation energy to the second band then a single band model is justified. Wave 
functions localized at an individual lattice site, w(x), which are called Wannier wave 
function are introduced and the energy eigenfunctions are expanded in the Wannier basis 



The Bose-Hubbard Hamiltonian follows from this expansion, where the hopping energy 
between matrix elements is given by 



In the last step the trapping potential is assumed to be approximately constant over the 
spatial variation of a single Wannier function. 

One expects a zero-temperature quantum phase transition from the superfluid (SF) 
state to the Mott insulator (MI) state and as the depth of the lattice is increased for 
integer fillings. The SF state supports long-range phase coherence while in the MI state, 
the atoms are localized and there is no phase coherence. This transition was observed in 
ultracold atoms by Greiner et al. (2002). 



\|/(x) =Y,biw(x-Xi). 




the on-site repulsion is 




and the energy offset due to the lattice is 
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1.6.3 Atom Chips 

Quasi-one-dimensional systems have also been realized on atom chips. The miniatur- 
ization and integration of matter-wave optics has led to the development of atom chips 
(Folman et ah, 2002; Fortagh and Zimmermann, 2007). It is now possible to confine, 
manipulate and measure atoms on a single device using electric, magnetic and optical 
fields. Bose-Einstein condenstation has been created in magnetic microtraps (Ott et ah, 
2001; Hansel et ah, 2001). The traps are highly elongated and the one-dimensional 
regime is realized when the transvere confining potential, fm^ is much greater than the 
relevant energy scales of the system, the thermal energy, ksT and chemical potential, 
/j. Esteve et al. (2006) realized both the ideal Bose Gas as well as the quasicondenstate 
in a quasi-one-dimensional trap on an atom chip. Theoretical work has investigated the 
transition from the ID Bose gas to the quasicondensate (Bouchoule et ah, 2007) as well 
as the growth of the quasicondensate (Proukakis et ah, 2006). Other experiments in 
one-dimensional traps include the demonstration of the first phase-preserving matter- 
wave beam-splitter on an atom chip (Schumm et ah, 2005) and of an atom Michelson 
interferometer on an atom chip (Wang et ah, 2005). 

1.6.4 Classical Field Model of Bose Gas 

In this thesis we numerically study the Bose-Hubbard model, presented earlier, within 
the classical-field approximation. The classical field approximation is equivalent to the 
first-order mean-field approximation. In this section we outline the validity of the clas- 
sical field approach for studying the dynamics of interacting Bose gases. 
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The dynamics of a BEC can be well-described by the Gross -Pi teavskii equation 
(GPE) (Gross, 1961; Pitaevskii, 1961) 

m *J!M_ = (_^ V + V(r) +^(r,0| 2 ) ¥(r,f), 

where the coupling constant is given by g = 4nh 2 a/m (a = scattering length, m = mass). 
The Gross -Pitaevskii equation is a mean-field approximation and is equivalent to the 
continuous nonlinear Schrodinger equation which is integrable. The GPE has been used 
extensively to describe the dynamics of the condensate in three-dimensional systems. 

Several studies have looked at the applicability of the mean-field or classical field 
description beyond the dynamics of the condensate. Kagan and Svistunov studied the 
evolution of an interacting Bose gas from a strongly non-equilibrium state towards con- 
denstation. They demonstrated that the classical-field description accurately describes 
a weakly interacting Bose gas in the absence of a condensate provided that the occu- 
pation numbers of the initially occupied state are much greater than unity (Kagan and 
Svistunov, 1997). Given this condition, the time evolution of a state can be accurately 
described by the diagonal elements of the statistical matrix in the coherent state repre- 
sentation. 

Castin studies the classical field model for one-dimensional weakly interacting Bose 
gases (Castin, 2004). The classical field model is generated by replacing the quantum 
mechanical operator \jlr(z) with a complex field V|/(z). For the interacting Bose gas, the 
state of the classical field is governed by a single parameter, 

^_ h 2 p 2 pg 
mksT kgT 

where p is the mean density, T is the temperature, m is the mass, and g is the interaction 
parameter. Castin calculates the correlation functions gi(z) = (\^(z)fy(0)) and g2(z) = 
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<\)/ t (z)^ t (0)\)/(0)^(z)). The contrast, C = g 2 (0)/g 2 1 (0) drops off quickly as % increases 
and then slowly approaches unity for % 3> 1, indicating that density fluctuations are 
suppressed for large %. 

The conditions for the validity of the classical field model for % 3> 1 are summarized 

as : 

1. large occupation numbers, k B T 3> /j 

2. gas is degenerate, pX 3> 1 

3. weakly interacting regime, 3> 1 

where £, = {h/m^i) 1 / 2 is the healing length. Note that the condition that x> 1 automat- 
ically satisfies conditions (2) and (3). 

In summary, the classical field model is a good approximation for weakly-interacting 
particles of a degenerate gas for large occupation number, in which case fluctuations are 
suppressed. 

Mishmash and Carr study the correspondence between the mean-field and the fully 
quantum BHM in the dynamics of atoms in ID optical lattices (Mishmash and Carr, 
2008). The mean-field BHM is equivalent to the discrete nonlinear Schrodinger equation 
(DNLS). They numerically investigate the analogs of dark soliton of DLNS in BHM and 
use the time-evolving block decimation algorithm (TEBD) developed by Vidal (Vidal, 
2004) to carry out the full quantum calculations. 

1.6.5 Chaos and Integrability in Quantum Systems 

Access to one-dimensional systems of ultracold atoms in optical lattices has led to real- 
ization of some known integrable models and observed effects of integrability in the 
dynamics of these models. We focus on the effects of integrability in bosonic systems in 
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optical lattices. The Lieb-Liniger model is a completely integrable quantum description 
of one-dimensional bosons with two-body 8-interactions (Lieb, 1963; Lieb and Liniger, 
1963). The Hamiltonian is 

^ = -Ett + 2c E 5 (^-^), (1-13) 

i=\ dx i (ij) 

where the interaction, governed by the parameter c is repulsive. The Lieb-Liniger model 
has been solved via Bethe Ansatz. In the limit of infinitely strong 8— repulsions, c — > °°, 
the hard-core bosons, also known as a Tonks-Girardeau gas, map to non-interacting 
fermions (Girardeau, 1960). While these models were proposed a half-century ago, the 
Tonks-Girardeau gas was only recently realized experimentally in Rb 87 atoms that were 
strongly confined in two directions in an optical lattice to create one-dimensional tubes 
(Paredes et ah, 2004). By applying a shallow lattice in the longitudinal direction, the 
effective mass and thus interaction strength were increased in order to reach the Tonks- 
Giradeau regime. 

Later experiments observed the effects of integrability on thermalization in a one- 
dimensional Bose gas. Kinoshita et ah (2006) demonstrate the first experimental evi- 
dence for the lack of thermalization in a many-body system with a large number of 
degrees of freedom for bosons in optical lattices. A gas of interacting bosons was pre- 
pared out-of-equilibrium by applying a laser pulse to a one-dimensional Bose-Einstein 
condensate in an optical lattice. For both strongly- and weakly-interacting bosons, the 
expanded momentum distribution retains the initial double peak structure. Even with the 
background harmonic potential, the system is integrable in the limit of infinite-strength 
repulsion. It was expected that a system with finite interactions, which is believed to 
be non-integrable in the presence of a harmonic trap, would reach thermal equilibrium. 
However, the absence of thermalization occurred even for finite interactions. Figure 1.5 
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Figure 1 .5: f(p e \), momentum distribution for (a) Yo=4, x=34 ms. tt>i ue = 15x, t re d — 30x 
(b) y =l, x=13 msJuue = !5x, t red = 40x (c) y =0.62, x=13 ms. t blue = 15%, t red = 40x. 
Green is the initial momentum distribution averaged over the first period. (Kinoshita 
et al, 2006) 

shows the expanded momentum distribution for three different coupling strengths. From 
the three peak structure in Fig. 1.5(a)-(b) it is clear that the gas has not thermalized for 
the Tonks-Girardeau limit (yo = 4) and the intermediate regime (Yo = 1) even after thou- 
sands of collisions have occurred between atoms. 

1.6.6 Constrained equilibrium 

One question that arises from these experiments on hard-core bosons is: do integrable 
systems, which don't relax to the usual thermodynamic equilibrium distribution attain 
some other steady state? Numerical studies on one-dimensional hard-core bosons in a 
lattice addressed the relaxation dynamics of a fully-integrable quantum system (Rigol 
et al, 2007). 
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The method to derive the steady-state distribution is to maximize the entropy subject 
to the constraints of the system, which include all of the conserved quantities (Jaynes, 
1957a,b). In this approach the many -body density matrix is given by 



p = Z 1 exp 



where {I m }, {'km} are the Lagrange multipliers which are determined by the initial con- 
ditions. This distribution is called the generalized Gibbs ensemble or fully-constrained 
thermodynamic ensemble. 

One-dimensional hard-core bosons on a lattice can be mapped to free-fermions via a 
Jordan-Wigner transformation. The conserved quantities are moments of the fermionic 
momentum distribution. Rigol et al. solve analytically for the density matrix with the 
constraints from the fermionic momentum distribution. The results of numerical simu- 
lations confirmed that when the system is prepared in the ground state of a small box and 
then allowed to expand in a larger box, it reaches a steady state, which is in agreement 
with the analytic results for the fully constrained system, rather than the grand canon- 
ical thermodynamic distribution. Additionally for an initial state with two momentum 
peaks, the two peaks structure remains after many oscillations, in agreement with the 
experiment performed by Kinoshita et al. (2006). 



1.7 Outline of Thesis 

In this work we present a comprehensive study of chaos and thermalization in the ID 
Bose-Hubbard model within the classical-field approximation. We study the threshold 
for chaos and its relation to thermalization. Two quantitative measures of thermaliz- 
ability are compared: the Finite-time Maximal Lyapunov exponents (FTMLE) and the 
normalized spectral entropy (NSE). The FTMLE, averaged over phase space, converges 
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to the maximal Lyapunov Exponent, the standard measure of chaos. A positive MLE 
indicates that points that are initially close in phase-space diverge exponentially, rather 
than linearly. The spectral entropy measures the distance between the time-averaged 
momentum distribution of the numerical results and the momentum distribution pre- 
dicted by thermodynamics, within the independent mode approximation. We investigate 
the dependence of the averaged FTMLE and normalized spectral entropy on a dimen- 
sionless nonlinearity parameter and the energy-per-particle, both of which are finite 
in the thermodynamic limit. The BHM is found to have a threshold for chaos which 
depends on the nonlinearity and the energy-per-particle. We study the size scaling of 
the Lyapunov exponent and normalized spectral entropy. 

Furthermore we study resonances in the Bose-Hubbard model to find the Chrikov 
criterion for chaos. The criterion predicted by the Chrikov criterion is different from the 
one inferred from numerical calculations, signifying the failure of the standard Chirkov's 
approach. 

There are at least three near-by integrable models: the Ablowitz-Ladik lattice, the 
continuous nonlinear Schrodinger equation and the noninteracting model. We outline 
the method of Inverse Scattering Transform and generate all of the integrals of motion 
of the closely related, fully integrable model of Ablowitz-Ladik. Furthermore, we dis- 
cuss the possible role of these conserved quantities in relaxation in the BHM. We con- 
jecture that the presense of quasi-conserved quantities may alter the scaling of the chaos 
criterion. 
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Chapter 2 

Thermalization and Chaos in the ID 
Bose-Hubbard Model 

2.1 Introduction 

One of the fundamental assertions of statistical mechanics is that the time average of a 
physical observable is equivalent to the average over phase-space, with microcanonical 
measure. A system for which this is true is said to be ergodic and one can calculate 
dynamical properties of the system from static phase-space averages. While this is 
believed to be true, because of the success of statistical mechanics in accurately pre- 
dicting experimental results, many open questions remain. Is ergodicity sufficient to 
ensure the accuracy of statstical mechanical predictions for times that are relevant for 
observations? 

2.2 BHM: Hamiltonian and Equations of Motion 

We study the dynamics of an interacting one-dimensional Bose gas on a lattice (ID 
Bose-Hubbard model (BHM)) (Jaksch et al, 1998) with periodic boundary conditions 
in the classical field approximation. The Hamiltonian of the system of interest can be 
studied in many different forms through canonical transformations. Each equivalent rep- 
resentation has a different Hamiltonian, canonical coordinates and equations of motion. 
A well-choosen canonical transformation can pose the problem in a way which is more 
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intuitive. Here we present three different representations: The real-space representa- 
tion, the momentum- space representation in terms of classical fields and the action-angle 
representation. 



2.2.1 Real-Space Hamiltonian 

In real space, the Hamiltonian is 



H = -J^s (¥.v + i + ^£l¥.s| 4 . (2.1) 

s L s 



The equations of motion are given by 



j t Vs = -J^: = ~ [ [- J (Vs+i +Vs-i -2%) +voN s \Vs\\s] (2.2) 



^f? = /[-7(v:+i+¥:-i-2v.:)+A'o^ s iv.vi 2 vj/.:], (2.3) 

and the canonical pairs are Qj = 2 S = ih\\r*. These equations of motion are equiva- 
lent to the discrete nonlinear Schrodinger equation (DNLS). The time-evolution of the 
fields is carried out in real-space in all of the calculations. 

2.2.2 Momentum-Space Representation 

Another set of canonical coordinates, the momentum-space fields, \\f n = \\f(k n ), are 
related to the real- space field \j7 s = ^f(x s ), by 



1 W S 

\\f n =—Y,Vse~ lkn ' Xs 

V ™s i— 1 
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where k n = 2%n/L and x s = sa. 

In the momentum- space representation, the Hamiltonian is 

// = Wto m | ¥m | 2 -^| ¥m | 4 ) +^£ ¥/ >; W/+; _ j (2.5) 

where the sum carries the restrictions: j ^ /,/. Several canonical transformations have 
been performed in order to write the Hamiltonian in this form. 

The canonical pairs are Qn = V|/ n , fp„ = z/U|/*. and the equations of motion are given 

by 



with the restrictions j ^ n,i on the sums and the indices span the range n, z, j = 
0, ±1, ±2, . . . , ( N * is supposed to be odd). 

The bare frequency of each momentum mode is given by 



=0— SOW 1 — &)>■ <2 - 8) 



The coupling constant is /jq = UNa/N^. Here J and U are the nearest-neighbor site- 
hopping and on-site repulsion constants of the standard Bose-Hubbard model respec- 
tively, and N a is the number of atoms. 

Time propagation is performed in real space, while the output and analysis of the 
numerical calculations focus on the momentum fields. 
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Table 2.1: List of Parameters and Variables 
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= Talbot time 


[/ 


= B-H on-site repulsion energy 


M>i 


= ground state energy of non-interacting model with quadratic dispersion 


5 


= size-dependent nonlinearity parameter 



2.2.3 Action-Angle Representation 

Equivalently, the momentum- space Hamiltonian can be written in terms of action-angle 
variables, by performing serveral canonical transformation on the momentum-space rep- 
resentation. It is this action-angle representation that will be the starting point of the 
resonant approximations and studies of individual resonances. The Hamiltonian is 

H = £ U m I m - + » £ (l mWj ) 1/2 e <«^-e,-e,) (2 . 9) 

L!l 1(1 m,l,i,j 

where the sum carries the restrictions: m + 1 = i + j; m ^ i, j; I ^ i, j. The momentum 
wavefunction canonical variables, {(v^, ih\f* k )} are related to the action-angle canonical 
variables { (4, Q k ) } by \|/^ = W \ e~' 6k ,\\f* k = J \e® k . In this form, the Hamiltonian can 
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be seen as the sum of an integrable term and a perturbation, H({Ik, 0^}) = //o({4}) + 
V ({4, where the integrable Hamiltonian is 



Ho = £ (to m 7 w - = £ (fao m |\|/ m | 2 - y IVml 4 ) • (2.10) 



Throughout the wavefunction \|/„ is normalized to unity: |\|/„| 2 = 1. 



Dimensionless nonlinearity parameter We define the dimensionless nonlinearity 
parameter, 

J J(N a /Ns) 



whose physical meaning is the ratio between the typical interaction energy per site 
U(N a /N s ) 2 and the kinetic energy per site JN a /N s . Note that this parameter governs 
both the strength of the nonlinearity and the strength of the perturbation from the inte- 
grable Hamiltonian (2.10). 

In Tables 2.1 and 2.2, variables are listed and the relationship between relevant 
parameters are summarized. 



2.2.4 Validity of the Classical-Field Theory 

Based on the studies of the validity of the classical-field theory for Bose gases (Castin, 
2004; Kagan and Svistunov, 1997) discussed earlier, the classical-field approximation 
will apply for the lattice site occupations satisfying 

— 3> max(K, 1) max 
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Table 2.2: Relationship between Parameters 
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M)i \ 27 V 



h 2 



1 — cos 



/ 2%n 

\n7 



27 



1 — cos 



(2%n 

\n7 



where An is the typical width of the momentum distribution. We note that the Mott 
regime, 7V a = integer x N s , An = N s , U/J > 2.2N a /N s (Hamer and Kogut, 1979), lies 
well outside of the above criteria. 



2.2.5 Nearby Integrable Models 

There are several known intergrable models that are limiting cases of the BHM. These 
include 
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1. Continuous Nonlinear Schrodinger Equation: In the contiuum limit, the ID 
BHM becomes 

h 2 . 1 „. r L . 1 



2. Linear Model: In the non-interacting limit, K — > 0, the ID BHM becomes a sum 
of harmonic oscillators 

= £fao m |\|/ m | 2 . 

m 

3. Independent Mode: If the interating term of vanishes, the ID BHM becomes a 
sum of decoupled nonlinear oscillators, 



H = £ftco m |y m | 2 - ^-\Wm\ 4 



with nonlinear frequencies given by Q. m = h(O m — ^|\|/ m | 2 



4. Ablowitz-Ladik Lattice: An alternate discretization of the NLS yields the 
Ablowitz-Ladik lattice, with Hamiltonian, 

^ = -E(^+i + ^+0 -^E ln ( 1_ ?l^l 2 ) ' 

n ° n V Z ' 

which will be discussed further in later chapters. 

2.3 Time Dynamics 

First we study the time dynamics of the ID BHM on a lattice with Af s = 21 modes. The 
system is prepared in a state that is narrowly distributed in momentum space and evolves 
according to the classical equations of motion. Initially the lowest three momentum 
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modes are occupied, the minimum number of modes required by selection rules for 
non-trivial processes leading to population of initially unoccupied modes. In Fig. 2.1, 




Figure 2.1: Time evolution of the wavefunction at the center of the box for a typical 
run for K = 0.009,0.09,0.36,0.9 with identical initial conditions, (a) Time dynamics 
of the real part of the wavefunction, Re \]/ JC= o (b) Frequencies of the real-space 
wavefunction, |FT[\j/ x= o(?)] | 2 . For K = 0.009,0.09 the descendants of the unperturbed 
frequencies generated by the first, "integrable"term of the Hamiltonian (2.5) are labeled. 

the time dynamics and power spectrum of the wavefunction at the center of the box, 
y x= o(t), are plotted for various interaction strengths for a typical initial state. As seen in 
Fig. 2.1(a), the time evolutions of the zero momentum mode is quasi-periodic for weak 
interactions, with a few easily identifiable frequencies entering the dynamics, which is 
confirmed by the power spectrum in Fig. 2.1(b). As the nonlinearity increases, more 
frequencies determine the dynamics and for sufficiently large nonlinearity the motion 
loses its quasi-periodic character and appears to be chaotic. 

The clear distinction between quasi-periodic and seemingly chaotic behavior of the 
time dynamics leads to the following questions: 

1 . Is the motion really chaotic? 

2. If it is, where is the chaos threshold as one increases the nonlinearity K? 

3. When chaotic, does the system reach thermal equilibrium? 
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In order to answers these questions, it is necessary to define appropriate measures of 
chaos and thermalization. 



2.4 Chaos: Calculating Lyapunov Exponents 

The standard signature of the chaotic nature of a region in phase-space is that the sep- 
aration between trajectories that are initially close grows exponentially with time, for 
typical trajectories, as captured by a positive maximal Lyapunov exponent (MLE). In 
regular regions the separation grows linearly (Chirikov, 1979), resulting in zero MLE. 
As we increase K in our system, we expect the phase space to change from being dom- 
inated by regular regions for small K to being dominated by chaotic regions for large 
K. In the present section, we use the MLEs to quantify this transition to chaos, which, 
as we will see in the subsequent section, coincides with a relatively broad change from 
unthermalizability to complete thermalizability. 

Consider two trajectories x(t) and x{t) with initial points jco and xq, respectively. 
The separation 8x(t) = x(t) —x(t) initially satisfies a linear differential equation, and 
the duration of this linear regime grows without bound as the initial separation Xq — Xq 
goes to zero. The finite-time maximal Lyapunov exponent (FTMLE) corresponding to 
the phase-space point Xo (Eckhardt and Yao, 1993; Voglis and Contopoulos, 1994) is 
given by 



The limit % n — > °° gives the MLE, ^(xq), but the FTMLE are themselves of intrinsic 
interest (Eckhardt and Yao, 1993; Voglis and Contopoulos, 1994; Contopoulos et al, 
1978; Contopoulos and Voglis, 1997). We chose a convenient quantum mechanical 



^fin(*o) = Jim —In 



1 , ||*(ffin)-*(ffin)ll 




1*0 —xo 



(2.12) 
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metric, ||jc — jc|| 2 = |\|/„ — V|/„| 2 . This metric becomes Euclidian under the canonical 
transformation 



Qn 

p„ 

||jc-jc|| 2 

where the sum runs from n = 



- (2^) 1 / 2 Re( ¥ J 
= (2^) 1/2 lm( ¥ J 

^I[(q:-q«) 2 +(p:-p,) 2 ] 

(iV s -l)/2ton=(iV s -l)/2. 



2.5 Thermalization: Calculating Spectral Entropy 

In order to measure thermalization, it is necessary to make thermodynamic predictions 
and compare the dynamics from the propagation of the equations of motion, with the 
thermodynamic state. In this section a method for calculating the thermodynamic state, 
within the Hartree-Fock approximation, is laid out. A full account is given in Appendix 
A. Additionally the spectral entropy is defined, which is a quantitiative measure of the 
difference between the time-averaged dynamical state and the expected thermodynamic 
state. 



2.5.1 Conserved Quantities 

A treatment of the thermodynamic state must take into account all conserved quantities. 
The known conserved quantities of the ID BHM are the energy and norm. Conserva- 
tion of norm, L^ =1 |y n | 2 , is associated with U(l) symmetry in the real/imaginary plane 
represented by the tranformation — > {y„e !6 ,\|/*e~' e }. Unlike the continuous 
NLSE, the total momentum is not conserved. 
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2.5.2 Hartree-Fock 

Within Hartree-Fock the form of the density distribution function is taken to be Gaussian 
and the thermal expectation value of the Grand Potental, 

(F) = {H)-T{S)-v{N a ), (2.13) 

is minimized, where N a is the norm. The density distribution function with two-body 
interactions has the form, 

ghf = -^exp (^ -ha n rf\\f„Y n ^J = -^exp ^-£a„/„^ . (2.14) 

We use the independent mode approximation so that in the second step, the off-diagonal 
elements are taken to be zero. The a n coefficients are unknown and are determined by 
the condition of minimizing the grand potential. The density distribution function is 
normalized, so that the integration of Ohf over all of phase space is 1, by 



Z = —L- f d N B f d N Ie-^ a "'« = FT-?-. (2.15) 

The expectation value of each term in the Grand Potential is calculated, using the 

Hartree-Fock density distribution function Ohf- The expectation value of a generic 
observable is given by 

(O) = / ^ e / d N IO({I n ,Q n })a HF (2.16) 

= — L^fjai f d N B [ d N IOe-^ anI " (2.17) 

The expectation values of the relevant observables are listed below. 
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Norm 



£ 1 



Hamiltonian - Kinetic Term 



N co„ 



n=l a " 



Hamiltonian - Interaction Term 



<*) = fe E (w^ 1/2 W^-' ,9 " +e - e - e 'A = SE^-^ ("O) 

\ z, ' t m,p,#,r / " m,p u "i U P 

Entropy 

5 = ~ (271ft)* / I dNl ° HF log ° HF = Ns ~ E lo §( na i) (2.21) 

2.5.3 Minimization of the Grand Potential 

The thermal expectation value of the Grand Potential within Hartree-Fock is given by 

(F) = E 5. + » E J_^_ _ T U _ £log(taJ ) _„ E ^_ (2 . 22) 

Taking the variation with respect to a n , and setting it equal to zero gives 

fl = .^. 2 ^ L ± + L + "o. (2 . 23) 

Using £ m a m ' = hN a and solving for a n , 



<*n=j^ [tUBn + 2/JoNa ~ fl] (2.24) 
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The thermal expectation values of the occupation of momentum mode n become 

(In) = — = 7 — • (2-25) 

CC„ h(D n + 2/J Na-ILl 

In general, the coefficients fi and T are unknown and are determined by imposing con- 
straints on the norm and energy, which come from the dynamical code. The constraints 
are 

(2.26) 

e t ={h) = £>„</„> + g £(/ m ) (/„) = +^N 2 a 

n 1 m,n n 

Beginning with the expression for (/„), we can solve for T in terms of n, N a and energy, 



T = ®n{In)+2^N a {I„) -^(I n ). (2.27) 



Summing over n, 

T = ^-[E k + 2/jqN 2 - vN a ] (2.28) 

where = L« This expression for T can be substituted back into the contraints 

to reduce the system to two equations with two unknowns. Using the expression for 
temperature and normalization condition, a single constraint remains to be solved, 

The Hartree-Fock approximation is known to overestimate the interaction energy in 
the regime of strong interactions. For sufficiently large /jq, the Hartree-Fock interaction 
energy, ^qA^ becomes greater than the total energy resulting in negative kinetic energy, 



46 



where the kinetic energy is E k = Ej —fioN^. For this reason, we determine the temper- 
ature T and the chemical potential fx using the time-averaged numerical kinetic energy 
(along with the norm) instead of the total energy. The quantity E k in the thermal distri- 
bution is fixed to the time-averaged kinetic energy of the final state from the dynamical 
code. We fix the norm to its numerical value, and subsequently solve iteratively for the 
norm to find all parameters. In this way, the total energy is never used in the constraints. 

Additionally the solutions must satisfy the physical constraint that /„ > for all n, 
which leads to bounds on fj. For T > 0, the condition such that the denominator is 
greater than zero for all n, is tua n + 2fioN a —/j>0, which leads to an upper bound for 
fi, fi < 2fioN a . There is a critical kinetic energy that corresponds to infinite temper- 
ature, which leads to equal population of all the modes, (/„) = hN a /N s . The critical 
kinetic energy, which separates the positive and negative temperature regime can be cal- 
culated as Ek-cr = E„ ^h® n = f E„ 27 [1 - cos (2™)] = 2JN a . For E k > E k _ cr the 
temerature is negative, and the lower bound on fi is h(D n + 2fioN a — n < for all n or 
/j > 4J + 2/JoN a . Close to the critical kinetic energy, both the temperature and the chem- 
ical potential diverge. By expanding the norm in powers of t0 M /(// — 2fjQN a ) , an estimate 
for the chemical potential when E k = E k _ cr ± £ is n m 2/jq + E k _ cr ± £f_ CT / (2e) . The 
temperature and the chemical potential were computed individually for each initial con- 
dition used. 

In Fig. 2.2, the initial and time-averaged momentum distributions of a representa- 
tive state are plotted for k = 0.09,0.36 and 0.9, along with the thermal Hartree-Fock 
predictions, (\y n \ 2 ) = (T/N a )/(h& n + 2fi N a - ju). 

2.5.4 Spectral Entropy 

For coupled anharmonic oscillators, as in the FPU study, energy equipartition among 
the normal momentum modes signified thermalization. In the BHM, the additional 
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Figure 2.2: Initial, final and Hartree-Fock thermal momentum distributions for K = 
0.09,0.45, 1.8, starting from the same initial state. N s = 21. The initial state is a repre- 
sentative state and the final state is time- averaged. Zj is the total energy per particle. 

conservation of the norm modifies the quantity that is equipartitioned. To determine 
the best measure for the equipartition we use the variational Hartree-Fock Hamilto- 
nian (Castin, 2001; Ohberg and Stenholm, 1997), H UF = £„ h(£>„ F \\\f n \ 2 , where the set 
of Hartree-Fock energies {h(£>„ F } was regarded as the variational field. This procedure 
gives h(i>„ F = h(£> n + 2/JoN a —iJ, where /j is the chemical potential. 

The new quantity to be equipartitioned is the distribution of the Hartree-Fock energy, 

| ¥ „(»|%C 

qn() E».|¥„'(')l 2 fi< 

A quantitative measure of the distance from thermodynamic equilibrium is the spec- 
tral entropy 

S{t) = -Y,q n {t)\nq n {t). (2.31) 
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In thermal equilibrium, q n is equipartioned, maximizing the spectral entropy at a value 
Smax = log(iVs). A more convenient quantity to study is the normalized spectral entropy 
(Livi etai, 1985), 

•q(f) = 5 max-^(Q ^ (232) 

which is unity at t = and vanishes as the system approaches thermal equilibrium. 



2.6 Results: N=21 Sites, Three-Mode Initial Conditions 

Initially, we study the FTMLE for a class of initial conditions where only the k = 0, ±1 
modes are occupied. In this subspace we sample uniformly from the intersection of 
the microcanonical shells in energy and norm; the energy is chosen to be the infinite- 
temperature energy of the subsystem, and the norm is 1. For each value of K, we sample 
100 points, which we set as the initial points Xq. To each initial point we add a small 
random vector, as little as machine precision allows, to obtain the corresponding JCo's. 
Each pair we propagate for a time ? fln , short enough to ensure linearity of the evolution of 
8x(t) but long enough to be able to clearly distinguish chaotic trajectories from regular 
ones on a plot of ln&t(?) versus t: the former are straight lines of positive slope, while 
the latter are logarithm-like (Contopoulos et al, 1978). We also verify that the average 
of the FTMLE's over the ensemble of initial conditions does not depend on ? fin as long 
as both criteria above are satisfied. In Fig. 2.3 the averaged FTMLEs are plotted as 
a function of the interaction strength. There is a distinct regime with zero Lyapunov 
exponent for small K < 0.5 and a strongly chaotic regime for K > 1 where all initial 
conditions have positive exponent. Next we consider the relation between chaos and 
thermalization in the system. 
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nonlinearity parameter, K 



Figure 2.3: Ensemble-averaged finite-time maximal Lyapunov exponent, X, and normal- 
ized spectral entropy, T|, as a function of the nonlinearity, K. N s = 21. 

In Fig. 2.3 the spectral entropy of the final time-averaged state, also averaged over 
100 initial states (drawn from the same ensemble that was used for the Lyapunov expo- 
nent calculation) is plotted for each value of K. For large nonlinearities, K > 1, the 
normalized spectral entropy goes to zero, indicating remarkable agreement between the 
final state and the thermal predictions. Note that this corresponds to chaos threshold 
observed previously. For K < .5 the normalized spectral entropy is above .5 signifying 
that during the time evolution the state of the system remains close to the initial state. 
In the Fig. 2.4, the normalized spectral entropy is plotted versus the FTMLE for each of 
the 100 individual runs for K = 0.36, 0.54, 0.72, 0.9. As seen in the plot, an individual 
initial state with larger FTMLE tends to have lower spectral entropy, i.e. to relax to 
a state which is closer to the thermal one. Beginning at K « 0.5, where the averaged 
FTMLE is substantially non-zero, some of the initial states thermalize completely. 
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Figure 2.4: Normalized spectral entropy of the final time-averaged state versus finite- 
time maximal Lyapunov exponent for each of the 100 initial condition used to compute 
the averaged value for K = 0.36, 0.54, 0.72, 0.9. N s = 21. 

Method of Calculating Spectral Entropy For an individual realization, we calcu- 
late /jq and T given the norm and the final kinetic energy. The spectral entropy is cal- 
culated for (1) the initial momentum distribution and (2) the final mometum distribu- 
tion, (\\\f n (t)\ 2 ), which is the time-averaged population of momentum mode n from the 
dynamical code. 

2.6.1 Fluctuations 

In order to confirm that the system is thermal when T] — > 0, we investigate the scaling of 
fluctuations of the kinetic energy for large K. 

For a system with N particles and volume V in thermal equilibrium, the fluctuations 
of extensive observables scale with the size of the system as V~ 1//2 , in the thermody- 
namic limit, N — > °°, V — > °°, N/V = const. Consider an extensive quantity, O. The 
temporal fluctuations of O are given by 



a 2 = (0 2 )-(0) 2 . 
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0(t), and the relative fluctuations are Co/ (O). Now consider 
two idential systems, A and B, with corresponding extensive observables A and B , 
such that O = A + B . Then 

a 2 = (O 2 ) - (O) 2 . 

= (0 2 A + 2 B + 20 A B ) - (0 A + B ) 2 

= (0 2 A ) + (0 2 B ) + 2 (0 A B ) - (0 A ) 2 - (0 B ) 2 - 2 (0 A ) (0 B ) 

= o 2 0a +o 2 0b + 2 (0 A B ) - 2 (O a ) (O b ) 

For a system in thermal equilibrium, the two parts of the system are decorrelated, 
so that (O a Ob) = (0 A ) (0 B ) and the last two terms cancel. For a general rescaling of 
the system N' = aN and the extensive observable (O)' = a (O), the fluctuations scale as 
o' = yfdoo an d the relative fluctuations scale as o' / (O) = a _1 / 2 Oo/ (O). In contrast, 
consider the case where two identical systems in identical initial states are concatenated 
in a regime where the behavior is regular. In this case, there will be strong correlations 
between the two parts of the system. In the extreme case of A (t) = B (t), then Oo = 
2<3o A = 2<3o B and the relative fluctuations o' / (O) will be constant, independent of the 
size of the system. 

For the system under consideration, the thermodynamic limit is taken by scaling the 
number of atoms and length as = CcA^, N' s = (xN s , while the interaction parameter 
jUo and the lattice spacing, a, remain constant (U, J = constant as well). In order to 
simulate the thermodynamic limit, the initial conditions are generated by concatenating 
a duplicates of the real-space wavefunction of the reference lattice. This is equivalent 
to generating an initial state with momentum modes \\f' k , given by \\f' Q = \\fo, \|/ ±ot = 
\|/±i, from the initial state in the N s q lattice, Generating the initial conditions in 
this way preserves the average energy per particle, in units of J = h 2 /2ma 2 . A small 
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Table 2.3: Thermodynamic Limit and Scaling 
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perturbation is added to the initial wavefunction to break the symmetry associated with 
the translational invariance. 

For the thermodynamic limit, we want to take the case where N a — > °°,L — > <*>, with 
N a /L = const, with the lattice spacing, a, and interaction strength, g, remaining constant. 
Consider the case where Ll = ocL, N' a — aN a . The scaling of relevant parameters is given 
in Table 2.3. 

We study the standard deviation of the fluctuations for systems with N s = 21 sites 
and oc=2,3,4. To compare fluctuations for different lengths, we calculate: 



g(N s ,N s0 ) 



o Ek (N s )/E k (N s ) 



a Ek (Nso)/E k (Nso) 



(2.33) 
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Figure 2.5: Reltative fluctuations in kinetic energy. Normalized standard deviation/mean 
for N s =2 1 , 42, 63, 84 lattice sites and Nq = 2 1 . Data points are for five sample runs with 
equivalent initial conditions for different lattice sizes and K = 2.69. 



where Oe(N s ) and Ek(N s ) are the standard deviation and time average of the kinetic 
energy for a chain with N s lattice sites. The reference lattice size is A^o = 21 and the cal- 
culation begins after the system has already thermalized. In Fig. 2.5 we plot o(N s ,N s q) 
as a function of a = N s /N s q for various lattice sizes on a log-log scale for K = 2.69. 

-1/2 

As can be seen clearly from the plot, the fluctuations scale as N s , indicating that the 
fluctuations are indeed thermal. 



2.7 Chaos Threshold for Different Lattice Sizes 



Let us start from the notion that the parameter K introduced in (2.11) is the only dimen- 
sionless combination of the parameters of the problem that remains finite in the ther- 
modynamic limit, N s — > °°, Na/N s = const, / = const, U = const. Curiously, the chaos 
threshold for N s = 21 is at K « .5, i.e. K ~ 1 . Another observation comes from a related 





(N s /N s0 ) 1/2 - 

^ — _ ^ 


k = 2.69 
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work (Villain and Lewenstein, 2000) on chaos threshold in NLSE with hard- wall bound- 
ary conditions. The authors find that the boundary between regular and chaotic motions 
of momentum mode, n, is given by (/Jo|Vn| 2 ) / (h(Oin) ~ 1, where h(D{ is the lowest exci- 
tation energy, e.g. the energy of the first excited mode in the case of the Hamiltonian 
(2.5). Assuming that the shape of the momentum distribution \\\r n | 2 as a function of n/N s 
should be fixed in the thermodynamic limit, the left-hand side of the above relationship 
also remains finite. These observations lead to a conjecture that the chaos criterion 
involves only intensive parameters and observables, i.e. those that are finite in the ther- 
modynamic limit. Our test for the above conjecture is based on the fact that for a chaotic 




Nonlinearity parameter, k 

Figure 2.6: Averaged Finite Time Lyapunov exponent, X[J], for three different system 
sizes, iV s = 11,21,41. For each K, the same energy-per-particle was used for each lattice 
size. The error bars represent one standard deviation. 

motion the majority of the trajectories cover the whole available phase space, and as 
a result the LE becomes (for a given set of parameters) a function of the energy only. 
This implies that for the same energy-per-particle and the same nonlinearity parameter, 
K, Lyapunov exponents for different lattice sizes should be similar. In Fig. 2.6 the time- 
averaged finite-time maximal Lypunov exponent is plotted for three different lattices, 
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A^ s = 11, 21, and 41. For each K, the same energy-per-particle (in units of J) is used for 
all three lattices. 

From the plot it is indeed evident that the LE is universal with respect to the size of 
the lattice and that the values of the LE for /V s = 1 1 already give a very good estimate of 
both the value of the LE and the threshold. 

2.8 Two Parametric Theory of the Chaos Threshold 

The universality observed above suggests the most relevant pair of variables for map- 
ping the chaos threshold, namely K and the total energy-per-particle, ej/J. In order 
to test these parameters, we independently vary the nonlinearity, K and the energy- 
per-particle, £t- For each K and £j, we generate ten initial states with microcanoni- 
cal weight in the reduced phase space of three (n = 0, ±1 ) or five (n = 0,±1,±2 ) 
momentum modes. It is necessary to generate initial states with five momentum modes 
n = 0,±1,±2 because there is an upper limit on the energy of an initial state with only 
three modes occupied. The finite-time maximal Lyapunov exponent and normalized 
spectral entropy are calculated for each of the ten realizations and averaged over this 
ensemble. For the rest of this work we will use the term Lyapunov exponent (LE) to 
denote the ensemble-averaged finite time maximal Lyapunov exponent and normalized 
spectral entropy (NSE) to denote the ensemble-averaged normalized spectral entropy, 
unless otherwise specified. The total data represents over 3200 runs. 

In Fig. 2.7(a) contour lines of the LE for N s = 1 1 are plotted versus the nonlinearity 
parameter and total energy-per-particle. One can observe an initial plateau in the LE for 
X < X c = 0.02, given by the solid line. The threshold depends on both the nonlinearity 
and the total energy-per-particle. Based on the parameter regime investigated, it appears 
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Figure 2.7: (a) Contour lines of the averaged FTMLE versus the nonlinearity, K, and 
energy-per-particle, £7- = (H — Ho)/N a , where H is the Hamiltonian (2.5), and Hq = 
—2J+ (1/2)^0 is the ground state value of H. The solid contour line corresponds to 
X c = 0.01. The diagonal solid line represents the set of energies and nonlinearities used 
in Fig. 2.6. (b) Contour lines of the averaged normalized spectral entropy versus the 
nonlinearity, K and energy-per-particle. Solid contour lines correspond to T] = 0.68 and 
T| = 0.36. For reference, the threshold line from the FTMLE in (a) is plotted (dashed 
line). N t = U. 

that the threshold persists for small K no matter how much energy is present. For small 
energies it is unclear whether the threshold will persist or vanish for K ^ 1 . 

After crossing the critical line the LE increases with uniform slope. The critical 
line resembles a hyperbola with the point of closest approach to the origin at (k, e T ) ~ 
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(0.5, 0.27), so that the hopping parameter J appears to be a relevant energy scale. This 
is probably not an accident: for £7- 3> J the dispersion law (D n begins to deviate from the 
(quadratic) dispersion law of the continuous NLSE with periodic boundary conditions, 
which is integrable. 

The normalized spectral entropy was calculated for the same set of data runs and is 
plotted in Fig. 2.7(b). There are two solid contour lines, at T| = 0.68 and r\ = 0.36. The 
second contour line at r\ = 0.36 follows closely the dotted line, which is the threshold 
from the Lyapunov exponent. It is apparent that the two plots have the same general 
features, and that there is a strong correspondence between the presense of chaos and 
thermalization in the BHM. 




Nonlinearity parameter, k 



Figure 2.8: Data points used for interpolation for contour plots in Fig. 2.7 superimposed 
on data for Lyapunov exponent. 

A few features deserve discussion. 

In the region of T] > 0.68, X = 0, which is enclosed by the first contour line in the 
NSE and the x- and y-axes, the system relaxes to a state that is closer to the thermal 
state even though the Lyapunov exponent is zero. The time dependence of the spectral 
entropy reveals that this relaxation takes place very quickly, after which the spectral 
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entropy remains flat for many fundamental cycles, indicating that no further relaxation 
will take place (See Fig. 2.11). This raises a few questions: To what state does the 
system relax in this region? It is possible to describe it by a constrained ensemble, where 
the constrained quantities are the conserved quantities of near-by intergrable systems? 
We will return to these questions later. 

For X 3> A, c , the region of strong chaos, the majority of initial states relax to the 
thermal state and all final states are close to the thermal. It is likely in this region, that 
full relaxation is not seen due to slow relaxation times and the spectral entropy would 
vanish for longer propagation times. 

We consider the two limits, K — > 0, £j ~ 7 and £r/7 — > 0, K m 3. In the limit of small 
K, the chaos threshold and NSE contour line r\ = 0.36 overlap and for K — > 0, £r > 0.67 
converge to a value that is independent of the total energy-per-particle. For the parameter 
region explored there is no indication that the threshold will vanish, even for very large 
energies. This suggests a dependence on the ratio of the nonlinear to linear terms, similar 
to the critical Reynolds number found in FPU (Livi et al, 1985). In the opposite limit of 
£t/J — > 0,K > 1.5 the behavior is quite different. While the Lypunov exponent is zero, 
there is significant relaxation in the momentum distribution. It is important to note that 
in the limit that £7- — > the initial state approaches the state where only the n = mode 
is populated, which is also the thermal state and thus the normalized spectral entropy is 
not well-defined in that limit. For this reason, data is plotted for e T > 0.027, the lowest 
energies simulated. However even for £j ~ 0.027, relaxation is visible in the momentum 
distributions. 

In Fig. 2.9 the standard deviation of the normalized spectral entropy is plotted along 
with the contour lines at r\ = 0.36 and r\ = 0.68 from Fig. 2.7(b). The contour line at 
T] = 0.36 follows closely the chaos threshold from Fig. 2.7(a). Far above the threshold, 
where the r\ — > 0, the standard deviation is also small indicating that most of the states 
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Nonlinearity parameter, k 



Figure 2.9: Contour lines of the standared deviation of the NSE versus the nonlinearity, 
K and energy-per-particle. Solid contour line corresponds to at T] = 0.36 and r\ = 0.68 
from Fig. 2.7(b). N t = 11. 

themalize, as expected. Below r\ = 0.68, (in the region bounded by the axes), relaxation 
is minimal and the standard deviation is small, indicating that most initial states will not 
thermalize. In the vacinity of the threshold (k = X c , which is close to r\ = 0.36), the 
standard deviation is larger. We conjecture that this is because there is a large spread in 
the amount of relaxation expected for different initial states with the same parameters 
and/or that some states have not fully relaxed due to insufficient propagation times as a 
result of multiple relaxation time scales. 

In Fig. 2.10 the normalized spectral entropy is plotted for three different lattice 
lengths, A^ s = 11,21,41 for the same energy-per-particle at each K, using the same energy 
values as in Fig. 2.6. While the main features are similar for the different lattice sizes, 
the size scaling of the NSE is not as universal as the scaling of Lyapuonv exponent. For 
small k's with the same energy-per-particle, more relaxation is seen in larger lattices. 
This suggests that the number of modes in involved in the dynamics may plays a role in 
relaxation. In addition the standard deviation of the NSE is larger than for the Lyapunov 
exponents. In contrast to the size scaling of the Lyapunov exponents, the variance of 
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Figure 2.10: Averaged normalized spectral entropy, T|, for three different system sizes, 
iV s = 11,21,41. For each K, the same energy-per-particle was used for each lattice size. 

the NSE increases with smaller chains. We conjecture that the large variance can be 
attributed to multiple relaxation scales. For example, for N s = 11,K = 1.5, individual 
states reveal that while most relax fully to the thermal state there is a single state that 
remains very far from thermal, which is the cause of the large variance. 

The large variance of the spectral entropy is the reason that there are more features 
in the contour plot of the normalized spectral entropy compared with the contour plot of 
the Lyapunov exponent. Repeating the simulations for longer times would likely smooth 
some of the features of the NSE contour plot, and decrease the variance in regions where 
it is currently large. 

2.8.1 Thermalization Times and Slow Relaxation 

The time dependence of the normalized spectral entropy is plotted in Fig. 2.1 1 for K = 
0.09 and £j = 0.0817, which is in the non-chaotic region. The initial, thermal and final 
momentum distributions are plotted in the inset. The time-dependent spectral entropy is 
calculated from a running average over the momentum distribution and plotted in units 
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of x t ai, the talbot time, which is the period associated with the lowest frequency of the 
non-interacting system with quadratic dispersion. After an initial relaxation during the 
first talbot time, the spectral entropy saturates and remains flat for close to 1000x ta i. The 
momentum distributions confirm that there is some relaxation, but that the state remains 
far from the thermal state. 




Figure 2. 1 1 : Sample time dependence of normalized spectral entropy, T|. K = 0.09, Zj = 
0.0817, N s = 21. Inset: Initial (dashed red), final (solid blue) and thermal (dotted black) 
momentum distributions. 

In Fig. 2.12 the time dependence of the normalized spectral entropy is plotted for 
two different initial states that both have K = 0.54 and total energy £j = 0.197. In 
Fig. 2.12(a) the normalized spectral entropy vanishes indicating that the state relaxes to 
the thermal state, which is also seen in the final momentum distribution. The normalized 
spectral entropy drops in several stages suggesting that there are multiple relaxations 
time scales. In Fig. 2.12(b), the state does not fully relax during the observed propa- 
gation time. After the initial relaxation, which is very similar to the previous case, the 
normalized spectral entropy slowly relaxes further but does not vanish in the observed 
time. Both states have a positive finite-time maximal Lyapunov exponent and thus are 
in the chaotic regime. These observations brings up several questions. Given states in 
the chaotic sea with the same total energy and strength of nonlinearity, why do some 
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Figure 2.12: Sample time dependence of normalized spectral entropy, T|, for two differ- 
ent initial states. K = 0.54, £ T = 0.197 for both. Af s = 21. Inset: Initial (dashed red), 
final (solid blue) and thermal (dotted black) momentum distributions. 

fully relax while other do not? Will these states fully thermalize for longer propagation 
times? What are the relevant time scales? What governs the slow relaxation times? 

Comparing the momentum distributions for both plots, the initial momentum distri- 
bution is almost symmetric in Fig. 2.12(a) so that the total quasi-momentum of the initial 
state is close to zero. Total quasi-momentum is not a conserved quantity of the BHM, 
although it is a conserved quantity of the noninteracting model, the continous model (in 
which case is becomes the true momentum) and the Ablowitz-Ladik discretization of the 
NLS. The total quasi-momentum is zero in the thermal state. For small k's, there is very 
little redistribution among the momentum modes and thus the total quasi-mementum 
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is well conserved. We call the total quasi-momentum a "quasi-conserved quantity" 
because it is not actually conserved in the BHM, but it conserved in the nearby inte- 
grate models and thus is expected to be conserved in the BHM when is "close" to one 
of the integrable limits. Proposed future work includes the investigation of the role of 
the conserved quantities of the nearby integrable models in the dynamics of the BHM. 

While these plots are sample runs, the pattern just described is observed in other indi- 
vidual runs for different values of K and e T in the chaotic region. Relaxation occurs on 
multiple time scales and the propagation times used in the simulations are long enough 
for the fast relaxation, but are not always long enough for the slow relaxation. For 
a given set of parameters, (k,£j) there are different slow relaxation times for different 
initial states. Insufficient propagation times are one possible reason for large variation of 
the individual NSE's observed in Fig. 2.9. In the strongly chaotic regime, it is expected 
that the the normalized spectral entropy will converge to zero for longer propagation 
times. However it is also possible that for X, « X c , (T| 0.36) some initial states will not 
fully relax, even for very long times. Furthermore, for 0.68 > T] > 0.36 it is likely that 
the variance will remain large. It is clear from Fig. 2.1 1 that some states do not relax, 
even for very long times. 

In summary, we have observed a threshold for chaos in the BHM, which depends on 
two parameters, the strength of the nonlinearity, K and the total energy-per-particle, 
£t/J- Far above the threshold, the state relaxes to the one predicted by statistical 
mechanics. Below the chaos threshold, we observe relaxation to a non-thermal steady- 
state. For small nonlinearities, k's the chaos threshold and absense of thermalization 
persist even for large energy-per-particle, e T ~ J. For regions just above the threshold, 
there are multiple relaxation times, with different intitial states relaxing on different time 
scales. These observations bring up several questions: What is the origin of the chaos 
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threshold? What governs the long relaxation times? Is the nonthermal steady- state 
affected by the conserved quantities of the nearby integrable systems? 
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Chapter 3 

Resonance Model and Failure of 
Chirikov's Criterion 

In this chapter we study a Hamiltonian coupling a single set of modes for initial states 
with few modes occupied. First we use a perturbation theory expansion that is valid 
for small nonlinearities. Second we study nonlinear resonances of the one-dimensional 
mean-field Bose-Hubbard model to predict the chaos criterion by Chirikov's method of 
overlapping resonances. 

3.1 Perturbation Theory Study of BHM 

To study a perturbation theory expansion we introduce the a prefactor £ in the driving 
term, which will be set to unity in the end. 

The full Hamiltonian in the momentum-space wavefunction representation is: 

// = £(^Q) m |y m | 2 -^|^ ra | 4 )+e^ £ Ym>YvVi'Vj' (3-1) 

m V Z / Z m ,j, :i , f 

where the sum carries the restrictions: m' + 1' = i' + /; m' ^ i' ', /; /' 7^ z -/ , /. We define 
the unperturbed Hamiltonian as 

#0 = £(faOm|Vm| 2 -y|Vm| 4 ) , (3.2) 
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which consists of decoupled nonlinear oscillators. The frequencies of the unperturbed 
Hamiltonian are 

Q n = C0„ - (/j /h)\\\f n \ 2 = (bin 2 - (/j /h)\\\f n \ z . 
The equations of motion of the full Hamiltonian are given by: 

^ n = -hW n = ~ l ^ n ~^ (3-3) 
= -ia n y n -ie^ V YiVWn+i-i (3.4) 

l,i;if=n,l 

Now we make the following assumption: (1) All modes except n are treated as 
independent oscillators with equations of motion given by the unperturbed Hamiltonian, 
\|// = — i£l\\fi which has solutions of = and (2) the nonlinear frequency of 

mode n associated with the unperturbed Hamiltonian is fixed, Q. n = const. With these 
restrictions, the equations of motion of mode n become 

|x|/ M = -/n n¥n -/£^ £ %^ J ^ n+I _ J e- l ^ + ^'- J -^ t . (3.5) 
3.1.1 Dynamics of an Initially Unpopulated Mode 

First we study the maximum value of initially unpopulated mode in the perturbation the- 
ory expansion. We consider the initial state where a block of modes from [—N ma x,N max \ 
are equally populated and study the time dynamics of mode Q = N max + 1 . In order to 
do this, we seek a solution of 

^-Vn + i&nVn=f(t) (3.6) 
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We want to find an integrating factor h(t), satisfying h{t) = i£l n h{t) such that 



^[\ Vn (t)h(t)] = h(t)^y n + ia n h(t)y n = f(t)h(t). (3.7) 
at at 



The solution is h(t) = exp(i£l n t). Let f(t) = T.aS( a ) exp(/A a ?) where 



E= E 

« i,rJ¥=n,i 



8 (a) = -^mjvn+i-j (3 - 8) 

A a = &l - £lj ; — Q„+/_y. 

Next we integrate over time 

£ dtj f (y n {t)e iant ) = ^£g(ay^ +A °> 
i|/ n (T) e ^ - v|/ n (0) = £ g( °° (^(^+Aa)x _ x \ (3 9) 

The time-averaged value of \\f n (t) is given by 
1 f T 



|\|/„| 2 = lim - / dt\Mf n {t)[ 



1 f T g(a)g(P) (V (A «- A e )f - e I '( A «- fi ») f - e -<( A p- n ») f - l) 
= lim — / —. r - 

°\ P (3.10) 
^ S 1 " ' " P - [5(A a - Ap) - 5(A a - On) - 5(A p - Q„) - 1] 



^ (Q„ + A a ) (Q„ + Ap) 



= V g(a) f (P) [5(A a - Ap) - 1] , 

^(a„+A a )(^+Ap) L 
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We introduce the dimesionless variables, 

2 N 



which is directly proportional to the nonlinearity parameter, K, but scales with the size 
of the system in the thermodynamic limit. We also introduce a new time scale, x = Qb\t. 

Occupation of mode N max + 1 for Quadratic Dispersion We are interested in the 
population of mode Q = N max + 1 when initially modes (—N max ,N max ) are equally occu- 
pied with population vfr. 

When only low momentum modes are excited, the true cosine dispersion laws can 
be approximated by a quadratic dispersion relation, (0„ = (h\n 2 and the unperturbed 
frequencies can be written as £li = (0/(7 2 — ^|\j/| 2 ). Dropping terms 0(Q in the denomi- 
nator, the equation of motion for \\fQ can be written as 

giiSlQ+Cll-Slj-SlQ+l-^t _ i 

Vfl (o = -e^.-^y I D +Q n — 

UI,|/|<G;#Z 2 + Z ~ j ~ 
e i(l 2 -f-(Q+l-j) 2 )* _ e -iQ 2 ^ 
¥e (T)=-£^ £ g 2 + / 2-2_ (g + / ,2 (3.11) 

e K2(i-Q)(j-Q)-Q 2 )^ _ e -iQ 2 x 

¥C(X)=_4? 2('-e)0-e) 

where we make the substitution / = j + i — Q'm the last expression. Consider the case 
where modes 0, ±1 are initially occupied with equal occupation, \j/ and Q = 2. The 
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nonzero terms in the sum correspond to modes (2,0 -H- 1, 1), (2, —1 -H- 1,0), (2, — 1 -h- 
0, 1). The time evolution becomes: 



\|/ 2 (t) =^e 



e i(n 2 +n -2ni)f _ ^ e ;(n 2 +n_i-n -£2i)f _ j 

h2 

^2 + ^0-2^1 + & 1-^0-^1 

! ' 4x - 1 \ „ _ -xf 1 + e~' 2t - 2e~ i4x 



(3.12) 



So that for N max = l,the time evolution of | V|/2(x) | 2 is 



/ m? *-? -6 /3 -cos(2x) -2cos(4x) 

|\i/ 2 (T)r = ^ z \|/° ' 



(3.13) 



Occupation of mode N max + 1 for Cosine Dispersion In the previous case, we 
assumed a quadratic dispersion, which corresponds to the free space. The true dis- 
persion of the BHM, which is a lattice model, is cosine, 



C0„ = -G>1^2 f 1 -N?COS 



/ 2nn 



(3.14) 



and the time evolution of mode population |\|/ 2 (x) | 2 becomes 



¥2(0 = 



^^co S (4^)A ( g ^, 2 (co S (4^,) + l-2cos(2^)) _ j 

SV I -cos(4n/N s ) - 1 +2cos(27i/AQ 



+ 2 



gifJVf (cos(4ji/JVi)-l) _ j 
-cos(47t/Af s ) + l 



Defining Ai 
l)/roi, the 



Ar 2 (cos(4rc/Ar,) + 1 - 2cos(27i/AQ)/a>i, A 2 = Af 2 (cos(47i/A^) 



|V2(T)| Z = ^V 



Ai 



+ 2- 



A 2 



Ai 



+ 2- 



A 2 



(3.15) 
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V2(T)| 2 = ^V 




(3.16) 



This expression for the time evolution of the modes is expected to be accurate for 
small values of Comparison of these predictions with numerics shows very good 
agreement for small 

3.1.2 Nonlinear Frequencies in the Real-Space Dynamics 

For small values of ^, the perturbation theory expansion accurately predicts the time- 
dynamics of the real-space wavef unctions. Next we consider the same model as the 
nonlinear coupling increases. In Fig. 3.1 we plot the modulus-squared of the Fourier 
transform of the real-space wavefunction at the center of the box, \|/ x=0 (f) , is plotted 
for N s = 21 and t, = 0.1, 1,2, 3 for an initial state with momentum modes \|/ w=0 ,\|/„ =±1 
occupied. Using a non-interacting model (no coupling between modes, but nonlinear 
terms associated with each mode), the predicted value of each of the resonances is Q. n = 
(i) n — ^I n . Taking /„ as the time-averaged value of I„(t) . The driving term is of the form 
e i(n m +n„-n ; )^ SQ additional frequencies are expected to be linear combinations of three 
other frequencies. 
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Figure 3.1: Frequencies of the real- space wavefunction at the center of the box, 
\FT[y x = (t)} 2 \ for £ = 0.1, 1,2,3. 



f2p = 2£1q — 

= f^o + ^l- — ^1+ 

Q§ = f2o — ^1+ — ^l- 
£ly = 2Q.{. — CIq 



3.2 Chirikov's Criterion 

Next, we study the Bose-Hubbard Model in the Born-Oppenheimer approximation 
(BOA). In the BOA, the populations of all but one mode are fixed. The fixed modes 
rotate in phase-space with constant frequency. The motion of the free mode is coupled 
to the the fixed modes and motion is governed by the equations of motion. 
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In the next section, we deduce the Chirikov chaos criterion by studying single res- 
onances within the Born Oppenheimer approximation. Within the resonant approxima- 
tion it is assumed that near a resonance, the resonance dominates the motion, so that 
driving terms can be studies independently. Within the BOA, we assume that the action 
variable of the single mode under study is small compared to the other modes, p,q,r, 
and that the other modes can be well-described by the integrable Hamiltonian, H$ that 
describes independent modes, 



Solving the equations of motion for the action variable, I m = — ^ gives I m = const = l h 



Thus the action variable of modes p,q, r is fixed and the angle rotates with constant 
frequency. 

Summary of assumptions: Throughout this section, we assume that 

1 . Single resonance approximation: we isolate individual driving terms of the Hamil- 
tonian and study the resonances of these models. 

2. Born-Oppenheimer approximation: The population of all modes, except the one 
under study are fixed. The fixed modes have frequencies £l m = (^d m — j . 

3. Quadratic Dispersion: For low-energy modes the dispersion (0„ = 
2a>i(Ay27t) 2 [l - cos(2%n/N s )] can be approximated by a quadratic disper- 
sion, a>„ = (bin 2 . 




while the equation of motion for the angle variable, 9, 



= T § ives 




(3.17) 
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4. Equal population of fixed modes: The populations of the fixed modes are taken to 
be equal in order to reduce the number of parameters. 

3.2.1 Classification of Resonances 

We classify the single resonances into three categories. We use the notation (p,q) — > 
(n, r), where the first two modes are the "feeding" modes and the second two modes are 
the "filling modes". Mode n is the mode of interest, whose dynamics are governed by 
the resonant Hamiltonian, while p,q and r are three different modes with fixed action 
variables. Additionally, the momentum indices must satisfy p + q = r + n. The reso- 
nances are classified according to the order of the exponent of mode n in the driving 
term. The three classes of Hamiltonians are: 

1. First-order Resonance (p,q) — > (r,n) 

The two feeding modes p,q fill two different modes, r,n. The first-order resonant 
Hamiltonian is 

H n = - + ( 7 » W") 1/2 cos (9„ + 9 r - 9 p - 9,) . (3.18) 

2. 1-R Resonance (p,p) — >■ (r,n) 

A single feeding mode p fills two different modes, r,n. The 1-R resonant Hamil- 
tonian is 

H n = ® n I n - ^I 2 n + W) 1/2 cos (9„ + Q r - 2d p ) . (3.19) 

The 1-R resonance differs from the first-order resonance only in the prefactor in 
front of the driving term. 
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Table 3.1: Resonance Parameters 



Order 


Modes 


Bare Detuning 


Mean Occupation 
of Fixed Modes 


First 


p,q^r,n 
p+q =r+n 


Aoi = p 2 + q 2 — r 2 — n 2 


h = (Wr? /2 


1R 


p,p^r,n 
2p = r + n 


Aoir = 2p 2 - r 2 - n 2 


I lr = (I p 2 I r )V2 


Second 


p,q — )■ n,n 
p + q = 2n 


A 02 = p 2 + q 2 -2n 2 





Order 


Drive Frequency 


Nonlinear Detuning 


First 


Vi = Q. p + Q. q -Q. r 

= (o p + (a q -(a r --y(I p +T p -Ir) 


Ai = (vi-co w )/wi 
= Aoi-^"i 


1R 


Vi r = 2£lp — £l r 

= 2to p -co r -^(2/ p -7 r ) 

LL 


Al r = (Vl r -©„)/©! 
= Aoir-^lr 


Second 


V 2 = Sip + Clq 

= (Op + <O q --y(I p +I q ) 


A 2 = (v 2 -2to n )/(2a 1 ) 
= Ao 2 /2-^/ 2 



3. Second-order Resonance (p,q) — > (n,n) 

Two different feeding modes, p,q fill the same mode, n. The second order reso- 
nant Hamiltonian is 

H n = iO n I n - ^I 2 + (l p I q ) 1/2 cos (20„ - Q p - Q q ) (3.20) 

In Table 3.1 the important parameters and definitions are listed for the three classes 
of resonances. The parameters include: the geometric mean of the fixed modes, 7, the 
bare detuning, Ao, and the frequency of the drive v. For the nonlinear detuning, the 
second equality holds when all the fixed modes have the same population. For a generic 
driving term of each type, the bare detuning, Aq is always negative for 1-R terms, always 
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Table 3.2: Bare Detuning: Generic Values 



Order 


Modes 


Bare Detuning 


First 


n + m,n + 1 — > n,n + 1 + m 


A m =(' M _|_ _|_ /„ i n2 

- (n + m + /) 2 -n 2 
= — 2ml 


1R 


n + m,n + m — >■ n,n + 2m 


Anw =2(n + m) 2 

— (n + 2m) 2 — n 2 
= -2m 2 


Second 


n + m,n — m — >■ n,n 


A02 =(n + m) 2 + (n — m) 2 
-2n 2 
=2m 2 



positive for 2nd order terms and can be positive or negative for 1 st order terms (see Table 
3.2). Next we will study the resonances in each class of Hamiltonians. 



3.2.2 First-Order Resonances in BOA 

Consider a first-order resonance where n + r = p + q; p,q ^ n,r. Fix the populations 
of the feeding modes, p,q and the filling mode r, so that the only variables are the 
population and angle of mode n. The frequencies of modes p,q,r are fixed to their 
values in the unperturbed Hamiltonian (3.2). Solving 0„ = ^ for the angle gives 

% = (co„ - t = €i n {I n )t. (3.21) 
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In this approximation, which we call the Born-Oppenheimer approximation (BOA), we 
keep the eight terms corresponding to (p,q) — > (r,n) from the full Hamiltonian. The 
Hamiltonian for first-order resonances in the Born-Oppenheimer approximation is: 



H n =CQ„/„ - ^/„ 2 + ^ (I n lpl q l r ) ^ COS {% + Q r - Q p - Q q ) 



2h 2 n h 2 

-iOjn ~ + (Wr) ^ COS (0„ - V? + ^ ) 



(3.22) 



where v = Q p + Q g — Now we let I\ = (I p I q I r ) 1 / 3 , divide by c&i, set ft = 1 and use 
^ = Aio/(M)i) to get 



/*„ = = - |/ 2 + 4^ (/! ) 1/2 cos (0„ - w + * i ) . 
(fll (Oi 2 



(3.23) 



We seek a canonical transformation to a rotating reference frame (7„,0„ — > /«,©«), 
where the new angle, 0„ is slowly varying. Introducing the type 2 generating function 



&=(Q„-Vt + fa)I n , 



(3.24) 



the new canonical variables are given by 



0n= ^ = 0„-V? + (|)i. 



(3.25) 
(3.26) 
(3.27) 



The Hamiltonian, which transforms according to h n = h„+^-, becomes 



fc„ = (co„ -Vi)/ml„ - |/^ + 4^ /2 /y 2 cos0„. 



(3.28) 
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The corresponding equations of motion are 



4 = -^ = 4^ /2 (/ 1 ) 3/2 sin9„ (3.29) 

dh 7 3/2 
e„= ^ = -A 1 -^/ n + 2^cos0 n , (3.30) 

" In 



where Ai = Vi — C0 M . 
Resonance Condition 

Resonance occurs at the stationary points, that is 



In 



9;; 



_ =0 4^(£") 1 / 2 if /2 sin§* = (3.31) 

-3 II 

ft*,=° ^ -A.-^ + a^cose-O. (3.32, 



The first condition is satisfied by (a) /* = or (b) sinG* = 0. For /„ = the phase is not 
well-defined and corresponds to a stationary point, even though (3.32) is not satisfied. 
The resonances of this model will correspond to taking sin 9* = and solving (3.32) for 

-Aj 1 n /2 -^ 3 n /2 = T2Ui /2 - (3-33) 
Squaring both sides and rearranging gives, 

Ai A? 



P n + 2j% + ^In - 4T 3 = 0. (3.34) 
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The solution to this cubic equation gives the fixed points, 



T* - 



I _ 

l n2 — 



f* 



2A 



3^ 

2Ai 1 



+A + B 



2 (A + B)+i^(A-B) 



(3.35) 
(3.36) 
(3.37) 



where 



A = 



B = 



1/2" 



1/2' 



1/3 



1/3 



(3.38) 



(3.39) 



The real and imaginary parts of the roots are plotted in Fig. 3.2 for I\ = 1/3 and Aoi = 
—4, which corresponds to (n + 2,n + 1) — > (n,n + 3), which is the negative detuning 
closest to zero. For the initial condition studied by the previous perturbation expansion, 
where a block of consecutive modes are initially occupied, the bare detuning is always 
negative. Thus, we still study the case of negative detuning although it is possible to 
have a first-order resonance with positive or negative values. The value of £ the separates 
regions with three real solutions to one real solution is given by (^^j + if = or 



Ai Aoi - %eh 



3^c 



3^c 



-h 



Aoi 
'2h 



(3.40) 



As will be shown below there is another relevant critical value of £, at lower values, so 
this critical value is defined as £ C 2. There are three physical fixed points (7* > for 
< ^c2 and one physical fixed point for £ > 



79 




Figure 3.2: Real and imaginary parts of the fixed points of first-order Hamiltonian 
(3.28). Aoi = -4, h = 1/3. 

Separatrix Width 

There are two separatricies: one is associated with h n = another passes through the 
fixed point (^* 3 , 0* ±7t) and is defined by the contour h n = h n (I* 3 , ±7l). 

In this section we calculate the maximum value of /„ for the separatrix defined by 
h n = 0. The maximum value of the separatrix occurs at 0„ = ±mn for integer m. We 
thus solve h n = for cos(0„) = ±1. 

-A 1 / ;i -i^ = T2^t /2 /,! /2 . (3.41) 

Squaring both sides and rearranging gives, 

§ + 4 L 5+4|/«-64i? = 0. (3.42) 
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Figure 3.3: Action-angle phase-space plots for variables /„, 0„ for the first-order reso- 
nant Hamiltonian (3.28). (n + 2,n + l) -> (n,n + 3) (A i = -4). 7i = 1/3. (a) £ = 1 (b) 
£, = 3.5 (c) £ = £ c i = 4.3 (d) £, = 7 > ^ C 2. The black contour line corresponds to = 
and the white contour line to h n (I n ,Q n ) = h n (I* 3 , —%). 

Solving this cubic equation gives the roots, 

4Ai 

Ise P l = -^r+A+B (3.43) 

4Ar 1 V3 
he P 2 = -^--(A + B)+iY(A-B) (3.44) 

4Ai L N y/3 , 
Ise P 3 = --^--{A+B)-i^-{A-B) (3.45) 
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where 



8(|)%4,i + w2/r((gW 1/2 



1/3 



B = 



3 \ 
3 



8 (|) +4/ 1-,6^(g) + 2H 



1 1/3 



(3.46) 



(3.47) 



The second critical value of £ is given by j + 2/j 5 = or 



Ai A i - ^ch 3/x 7 

sir = — — = 



5Cl 



Aoi 



(1 + 3^2)/! 



(3.48) 



The maximum value of the separatrix occurs at 0„ = ±71, (cos(0„) = — 1) for ^ < 
and at Q n = 0±27t, (cos(9„) = +1) for £, > L, c \. In Fig. 3.4 the real and imaginary parts 
of these solutions are plotted. 

From the fixed points of the Hamiltonian, solutions of h n = and sample contour 
plots, the phase-space of the first-order Hamiltonian can be characterized as follows: 

£, < £ c i There are two resonances. One is at (/„, Q n =) = (7* 2 , ±7t) with the separatrix 
defined by the contour at h n = 0. This resonance can be seen in the phase-space 
diagrams of Fig. 3.3(a)-(b). The height of the separatrix is given by I sep 2 which 
increases from zero at ^ = to a maximum value /j^-max = 2\/2I\ at ^ = ^ c 2 
as seen in Fig. 3.4. For I\ = 1/3 the maximum value is /,,2-max ~ 0.84. The 
second resonance is given by , 9„ = 0) . The separatrix is defined by the contour 
that passes through the saddle point (/* 3 , ±7t). This resonance diverges as £, — > 
and remains above 1 in this parameter range is thus not physically accessible. 
However, the lower separatrix comes into the accessible phase space as seen in 
Fig. 3.3(b).Due to the normalization, the occupation of /„ is bounded by one. 
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10 



sepl 
— Isep2 ■ 
-^sep3 

15 20 - 



"0 5 10 15 20 

Figure 3.4: Real and imaginary parts of the solutions of h n = for the first-order Hamil- 
tonian (3.28). A i = -4 (n+ l,n + 2) -»■ (n,n + 3), 7i = 1/3. For ^ < £ c i, / ?ep2 gives the 
separatrix for the the resonance at (/* 2 , 0„ = ±7c). For £ > £ c i, 7^ ep i gives the separatrix 
for the the resonance at , 0„ = 0) . 

Larger values of I n are plotted in Fig. 3.3 to gain a of deeper understanding of the 
phase space. 

£, = ^ci At the first critical point, the separatricies of the two resonances overlap as show 
in Fig. 3.3(c). At this critical point, the lower bound of the separatrix defined by 
h n = h n (I*^—%) touches zero. Starting with an arbitrarily small occupation, the 
entire physical phase space becomes accessible to I n once the two resonances 
touch. 

£ c l < t, < t, C 2 The two resonances at (/* l7 0) and (I* 2 , ±7t) remain, with the separatrix 
of the first defined by h n = and the separatrix of the second resonance defined 
by h n = h n (I* 3 , — 7t). The association of the separatricies have switched from the 
case where \ < ^ c i- 
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£ > ^c2 At the second critical point the fixed points (7* 2 , ±71 ) and (/* 3 , ±7t) merge and 
annihilate and a single resonance at (/„,0„) = (^,0) persists. Above there is 
a single resonance at (/ M ,0 M ) = (/* 1; 0), as show in Fig. 3.3(d) and the separatrix 
is defined by the contour h„ = 0. The height of the separatrix is given by T sep i in 
Fig. 3.4. 

3.2.3 1-R Resonances in BOA. CL p = CL q 

Consider a first-order resonance where the feeding modes are the same, n + r = 2p. In 
this case, we keep four terms from the full Hamiltonian, and the Hamiltonian for the 
first-order resonances (1-R) is 

H n = t0 n 7„ - ^l 2 n + ^ (l n I 2 p I r ) 1/2 cos (9 n + Q r - 2Q p ) (3.49) 

After making the transformation to a rotating reference frame, introducing V\ r = 2£l p — 
Q. r , fixing the populations of modes I p J r , and dividing by Q) M , the Hamiltonian, h n = 
H n /(d n becomes 



~h n = -Ai r I n - |g + 2^ /2 /f r /2 COS0„ 



(3.50) 



where Ai r = Vi r — Q) n . The corresponding equations of motion are 



/„ = -^ = 2^ 1/2 /t /2 sin0„ (3.51) 
d0„ 



3/L 



'-3/2\ 1/2 

^1 „ \ 



The Hamiltonian for the 1R resonances differs from the first-order Hamiltonian only 
by prefactors. The resonances of the 1-R Hamiltonian will thus be very similar in form 
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to those of the first-order resonances, and can be determined from the previous analysis 
by a proper rescaling of I\ . 



3.2.4 Second Order Resonances in BOA 

For the second order case, we consider a single filling mode n and restrict the resonance 
to modes that satisfy p + q = 2n. The Hamiltonian for a single second order resonance 
in the Born Oppenheimer approximation is 



H n = - + J^h {iplq) ^ COS (20„ -Q p - Q q ) 

= UnIn ~ J^In + jph ^ COS (29 n - V 2 t) 

where in the second line we set I p and I q to their unperturbed values and 



(3.53) 



v 2 = Clpilp) + Cl q (I q ) =& p + & q -^ {Ip+I q ) ■ 



Next we make a canonical transformation to a rotating reference frame (7 M , W — > /„, 0„), 
through the type 2 generating function 

<D = I(20„-v 2 * + (|>i)/„. (3.54) 

The new canonical variables are determined by 

d<Z> - 
d0„ 

3$ 1 

« = 3T = x(20„-v 2 ? + (^i), (3.56) 
dl n 2 
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and the Hamiltonian transforms according to H n (I n B n ) = H n (I, B n ) + ^ becomes 



H H =[a>n- ^v 2 ) /„ - J±$ + (7,4) 1/2 cos 20„ 



(3.57) 



Next we divide by &\, which is equivalent to rescaling time by a factor x = (bit, intro- 
duce £ = ji/o/M>i and set H = 1. Furthermore, we define 7 2 = (I p I q ) 1 / 2 as the geometric 
mean of the filling modes I p J q and A 2 = (^V 2 — Cfl„)/cQi. 



(3.58) 




The corresponding equations of motion are 



3/z 

/„ = --^ =4^/2 sin 20 n 

9„ = 3^ = - A 2 - + 2^7 2 cos 29„ . 



(3.59) 
(3.60) 



If we furthermore assume that the filling modes have equal populations, I p = I q = 1% and 
that the dispersion is quadratic, (0„ = n 2 (b\ then the detuning reduces to 



Ai = Wi (v " 2c0m) = i t (p2 + ^ 2 _ 2 " 2)&i _ + 7 ^ (3 - 61) 

= ^-^2. (3.62) 
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Resonances Conditions 

Unlike the first-order resonances there is a simple expression for the fixed points of the 
second order resonances. The conditions for the fixed points are /„ | ^ g*j = 0, Q n | ^ = 
0. 



0« 



4^/ 2 sin20: = O 



-A 2 -^ + 2£/ 2 cos20* = O 



(3.63) 



(3.64) 



1. Casel: /* =0. 

The condition for /„ = is satisfied by /* = and the condition on 0* such that 
0„ = is 



cos 20! 



A 2 A 20 1 



2^/ 2 4^7 2 2 

Note that if /„ = 0, the phase of the mode n is not well-defined. In this case, the 
motion of the action variable is always stationary. Thus a second-order resonant 
Hamiltonian can never populate an initially unoccupied mode. However, once a 
small seed is present, the mode can grow and enter the dynamics. 

2. Case 2: sin 20* = (cos 20* = ±1) 

The condition for /„ = is satisfied by sin 20* = and the condition on /* such 
that 0„ = is 



/* = -j ± 2/ 2 = +I 2 ±2/ 2 (for cos20: = ±1). 
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Classification of Stationary Points 

Following Tabor (1989) we outline a method for determining the stability of the sta- 
tionary points by linearizing about the fixed points. Consider a general second order 
differential equation, written as a pair of first-order differential equations, 

x = f(x,y) (3.65) 
y = g(x,y) (3.66) 

that has stationary points {(x* ,y*)\f(x* ,y*) = 0, g(x*,y*) = 0}. The stability of these 
points can be determined by linearizing about the stationary points. 




5*1* >y dy\ x ,y 




(3.67) 



The solution to these equations of motion, x = Ax, is x = c\e^ lt y\ + C2e^ 2f V2 where 
Xj is an eigenvalue of A with corresponding eigenvector vj. The stability of the fixed 
points is determined from the eigenvalues, Xj = aj + ibj, aj,bj G R. The fixed points 
are classified as follows, 

1. center: = 

2. spiral: ^ 0. Stable for < 0. Unstable for > 0. 

3. node: = 0. Stable for < 0. Unstable for a\2 > 

4. saddle point: b\,2 = 0, «i < 0, «2 > 

The linearized equations of motion for the second order Hamiltonian (3.58) are 
-ikkti -§t%,e* \f 4^ 2 sin20: 8^ M 7 2 cos20* 

%\m ikkfiJ V "5 -4^sin20: 



(3.68) 



88 



For each set of fixed points, we calculate the eigenvalues of the matrix, given by 
det | A — AI| = to determine the type of fixed point. 

1. Fixed Point 1: (7*,0*) = (O^cos" 1 (A 2 /(2$7 2 ))) 
The eigenvalues are given by 



I = ±4^/2 



a 2 Y 



1/2 



= ±2\ 



(2/ 2 ) 2 -(^ 



1/2 



(3.69) 



These eigenvalues are real for |A 2 | < 2^/ 2 . Given that the bare detuning, A02, is 
always positive for the second-order Hamiltonian (see Table 3.2) and the typical 
population of the fixed modes, 7 2 is also positive, this condition leads to a critical 
value of the nonlinearity parameter, 



A 02 
6/ 2 



(3.70) 



For ^ > ^ c this fixed point is a saddle point and the separatrix passes through this 
point. 

2. Fixed Point 2: (/*,§*) = (-A 2 /£ + 2/ 2 , mn),m = integer 



I = ±2v / 2^/ 2 1/2 



A 2 



-27 2 



1/2 



(3.71) 



For I* = — A 2 /£, + 2/ > the fixed point is a center, and for /* < it is a saddle 
point. We exclude the unphysical values /* < 0. The stationary points for /„ > is 
at cos(20 n ) — +1, are the physical resonances for the second order driving terms, 
with resonant value given by 



^res — 



A02 

2^ 



+ 3/ 2 . 



(3.72) 
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Furthermore the condition for the existence of this type of resonance is £ > £ c . 

3. Fixed Point 3: (/*,§*) = (-A 2 /£ + 2/, (2m + l)n/2), m = integer 

For cos(20*) = — 1, the value of /„ is always negative and thus not relevant to the 
current analysis. 

Separatrix Height 

The separatrix passes through the fixed point at /* = 0,0* = ^cos -1 (A 2 / (2^/2)) when 
0* is real. The contour plots of h n confirm that the separatrix passes through this fixed 
point. The maximum value of /„ along the separatrix can be found by solving for l n 
when h n = 0. Given 

h n = -A 2 I„ - |S + 2^4/2 cos 20„ = 0, 

I n = - ^ + 47 2 cos 20 M = - ^ + 27 2 ( 1 + 2 cos 20 n ) . 



The maximum height of the separatrix occurs at 0„ = and is given by 



A02 



'sep 



6/ 2 . 



(3.73) 



In Fig. 3.5(a)-(b) the phase plots for the second order Hamiltonian is plotted for 
nonlinearities above and below the critical value. For ^ < there is no resonance in 
the phase-space that corresponds to physical values of l n . In Fig. 3.5(b), where ^ > t, c , 
there is a resonance and the separatrix is defined by the contour h n = which is plotted 
in black. The resonance values and separatrix height are labeled by the values given by 
(3.72) and (3.73), respectively. 
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Figure 3.5: Phase-space plots for the second order Hamiltonian. (a) £, = 0.5 (k = 
0.16 for N s = 11). There is no resonance in the range < /„ < 1 for £, < £, c = 1. (b) 
£ = 1.5 ( K = 0.49 for AT* = 11). h = 1/3, A 02 = 2. For £ > £, c the resonant value is 
given by (3.72) and the separatrix with is given by (3.73). 

The values of the action variable at resonance and the separatrix height are plotted 
in Fig. 3.6 for the two second-order Hamiltonians as a function of the nonlinearity, 
\ for the two lowest bare detunings. The two detunings are A02 = 2, corresponding 
to (n+ l,n — 1) (n,n) and ^ c = 2/2/3, and A02 = 8 for (n + 2,n — 2) — >■ (n,n) and 
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Figure 3.6: Resonant values and separatrix width for second-order resonance for A02 = 2 
[(«+ 1, n- 1) ->■ (n,n)] and for A 2 = 8 [(n + 2,n- 2) (n,n)]. 7 = 0.33. Vertical 
lines: A20 = 2: £ c i = I/3/2 = 1, A02 = 8: ^2 = 4/3/2 = 4. The upper x-axis gives the 
corresponding values of K for N s = 11. 

£ c = 4/2/3. It is clear from the plot that for the second order resonant Hamiltonian there 
is a critical value of the nonlinearity such that below that value there are no physical 
resonances present. 



3.3 Failure of Chirikov 's Criterion 

The Chirikov criterion for the onset of global chaos is governed by the ratio of the width 
of the separatrix and the distance between the resonances for individual resonances. 
When the width of the separatricies becomes comparable to the distance between the 
resonances, the Chirikov parameter Ki satisfies the inequality, 

separatrix width AI 

K = = = — > 1, (3.74) 

distance between resonances / res i — / re s2 ~ 
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and the system is predicted to be chaotic. 

We should note that the study of resonances in the BHM is quite different from the 
case of the kicked rotor. In particular, it is not possible to independently vary the fre- 
quency of the drive and the strength of the drive in the BHM. In addition, in the canonical 
transformation to the rotating reference frame the action variable is unchanged, /„ = /„, 
so that the resonant values of /„ are not well- spaced in the original action coordinates. 
Furthermore in the BHM, for the first order resonances, a single driving term has reso- 
nances at multiple action values. 

Note that for a generic driving term, (n + p,n + m — p) — >■ (n, n + m), the bare detun- 
ing is given by 

Ao = (ft + p) 2 + (ft + m — p) 2 — n 2 — (n + m) 2 = 

= n 2 + 2np + p 2 + n 2 + m 2 + p 2 — 2np — Imp + 2nm 

(3.75) 

— [ft 2 + ft 2 + 2nm + m 2 ] 
= 2p 2 — 2mp, 

which is always constant. 

We present several ways to deduce the criterion governing the threshold for the 
mean-field Bose-Hubbard model and see that the various methods are in agreement. 

Onset of Second- Order Resonances First we consider the case of the second-order 
resonance, where the analytic expression is simple. 

From the previous analysis for the second order resonance, the separatrix width is 
given by, 

/se P = := ^ + 6/. (3.76) 
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The distance between resonances is proportional to £, 1 , 



A/ r e S = /res(Afc) -7res(A a ) = 



Ap a - Apfc 



(3.77) 



For the lowest two resonances, a naive comparison of the separatrix height and distance 
between the resonances gives, 



1 sep 



A/ n 



-2/5 + 6/ 
6/5 



(3.78) 



and chaos exists for 



(3.79) 



However, but looking at Fig. 3.6 , we see that once the second resonance appears, the 
two resonances overlap. Thus we take the appearance of the second resonance as the 
criterion for the second-order case, which gives 



c2 



31 



or 



(3.80) 



Overlap of resonances within first-order resonant Hamiltonian For the first-order 
resonance, the resonances of 7* 2 grow from zero and exist for any 5 > 0. However, 
the separatrix of the resonance is small and thus even though the resonances may over- 
lap, the population is still expected to be confined in a narrow region of phase-space. 
However, there are additional resonances at larger values of /„, which are initially inac- 
cessible and move down from above. For the first-order resonance we use the criterion 
that the two separatricies of the same first-order Hamiltonian overlap - i.e. all of phase 
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space is accessible to a single mode. The condition for being able to explore the entire 
phase-space for a given /„ is thus 

(l + 3^2)7i 



which is equivalent to 



/£> A ° 1 ,_ . (3.82) 

S " (1 + 3^2) 



For the lowest bare detuning, Aoi = —4, the criterion is 

(3.83) 



7£ > « 0.84 

b " (1 + 3^2) 



which is the comparable to the results for the second order resonance. 

Dimensional Analysis An alternate way of coming to this conclusion is to consider 
dimensional grounds. The starting point is to assume: 

1. Typical mode occupation: 7~ ^ 

2. Resonant approximation: only include resonant terms in equations of motion. 

3. Quadratic dispersion instead of cosine: G)„ = (h\n 2 

These are the same ingredients for the previous analytic analysis. A quadratic dispersion 
leads to translational invariance. The momentum distribution can be shifted with no 
effect. The resonant equations are invariant under a shift in 'n', thus the frequency scale 
is set by (b\. As a consequence of the resonant approximation, the equations of motion 
only couple neighboring modes. Thus, /uo and An can only enter as Thus the only 
parameters are no /An and c&i. The only dimensionless combination of these quantities 
is hq/ (h(h\An) = £/An, so this parameter must be what governs the Chirikov threshold. 
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Signficantly, all of these methods give the same functional form of the criterion and are 
identical up to numerical factors. 

3.3.1 Thermodynamic Limit 

Next we consider how %, scales in the thermodynamic limit, (N s — > °° while U = 
const, J = const, N a /N s = const). Additionally, the width of the momentum distribu- 
tion remains fixed, An/N s = const. In this limit, the nonlinearity parameter K is constant 
while £ ~ Ng. The typical occupation of the fixed modes is given by / ~ (An) -1 ~ 
(Ng) -1 We find that the Chirikov parameter, 

_ N 2 

Stni~I~j±~N s (3.84) 

diverges in the thermodynamic limit, predicting that there is no threshold in that limit 
and Ki is always greater than one indicating that the system is always chaotic. However, 
this is not what we have observed numerically. At a minimum Chirikov predicts the 
threshold to scales linearly with the size of the system and our numerics indicate that 
the threshold depends on parameters that are independent of the system size. 

3.3.2 Continuous Limit 

To check this threshold, we also consider the continuum limit, in which the length of the 
system, normalization and interaction parameter are fixed, while the distance between 
lattice sites goes to zero: L,N a ,/iio = const, a — > 0,N S — > °°. 
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In the continuum limit, the number of modes in the initial momentum distribution, An, 
must diverge, so that N a = An|\j/„| 2 remains constant. Alternatively we can say that 
An/N s is fixed. From either perspective, it is clear that %, vanishes in the continuum 
limit, predicting regular motion, as expected due to the known integrability of the con- 
tinuous NLS equation. 

This leads to the question: Why does the Chirikov analysis fail in this case? Several 
possible explanations are: 

1 . non-quadraticity of the spectrum 

2. break down of the resonant approximation 

3. interferences between resonances. 

First there are possible corrections due to the lattice and that we have assumed a 
quadratic dispersion instead of a cosine dispersion. We conjecture that the Chirikov 
analysis is incorrect because the the resonant approximation is wrong and off -resonant 
terms are significant. That is the motion is not dominated by single driving terms. It is 
also possible that the single resonance approximation is incorrect and a multi-resonance 
model is necessary for recovering the correct scaling. 
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Chapter 4 

Conserved Quantities of the 
Ablowitz-Ladik Lattice 

From the numerical work on the Bose-Hubbard Model, we have seen regions where 
although individual trajectories are chaotic the system does not relax to the expected 
thermal state. Additionally there are regions in the parameter space of nonlinearity, K, 
and energy-per-particle, £j, where the system relaxes even though chaos is not present. 
Additionally it has been found that the slow relaxation times for individual states with 
the same nonlinearity and energy-per-particle vary widely. What is the origin of these 
phenomena? 

There are at least three nearby integrable systems of the BHM. These are the 
noninteracting case (k = 0), the continuum limit, which is the continuous nonlin- 
ear Schrodinger equation and the integrable discrete nonlinear Schrodinger equation 
(IDNLS) also know as the Ablowitz-Ladik (AL) lattice, which we discuss in this sec- 
tion. Looking at the relaxation of individual realizations, it appears that initial states 
which have a higher quasi-momentum also have a higher final spectral entropy. The 
total quasi-momentum is not a conserved quantity of the BHM. It is however, a con- 
served quantity of all three nearby integrable systems. In some cases the normalized 
spectral entropy many not vanish because the system has a very slow relaxation time. 
In other cases it may be that it will never fully relax because of the nearby conserved 
quantities. Are the slow relaxation times governed by the conserved quantities of the 
nearby integrable systems? In the regions where it is clear that the steady-state is not the 
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thermal state, what governs the steady-state? This evidence suggests that the integrals 
of motion of the nearby conserved quantities play a role in the relaxation dynamics of 
the BHM. In this section, we derive the conserved quantities of the AL. 



4. 1 Integrable Discrete Nonlinear Schr odinger (IDNLS) 
Equation 



The continuous Nonlinear Schrodinger (NLS) Equation is given by, 



iq + q xx + G\q\ q = 



(4.1) 



where q = dq/dt and q x = dq/dx, is know to be completely integrable. The correspond- 
ing Hamiltonian is 



H(q,q*) = - f 1 
JO 



1 



*\2 



dx, 



(4.2) 



with canonical pairs q,q*. Periodic boundary conditions are assumed. The mean-field 
Bose-Hubbard model is a discretization of the NLS equation. In real space, the Hamil- 
tonian can be written as: 



H = -Y,(l n 1*n+i + 1nQn+l-2\q„\ Z ) + ^£ \q n \ 



(4.3) 



with the Poisson brackets 



{<lm,q*n} = i5m,n, {q m ,q n } = {q* m ,ql} = 0. 



(4.4) 
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The equations of motion, which are given by q n = {//, q n }, are 

iq n = -(q n +i +1n-l -2q n ) + o\q n \ 2 q n . (4.5) 

This equation, also know as the Diagonal Discrete Nonlinear Schrodinger (DDNLS) 
Equation, does not preserve the integrability of the continuous case. An alternate dis- 
cretization, which is integrable (Ablowitz and Ladik, 1976), has equations of motion: 

iq n = -{q n +l +q n -l-2q n ) + o\q n \ 2 (q n+ i +q n -i)/2 (4.6) 

and is suitably called the Integrable Discrete Nonlinear Schrodinger (IDNLS) Equation, 
also know as the Ablowitz-Ladik lattice (AL). This equation can be derived from the 
Hamiltonian (Scharf and Bishop, 1991; Herbst et al, 1994) 

H = -£ (q n q* n+ i +q*„qn+l) ~ ^£ln (l -% n \ 2 ) (4.7) 
with the nonstandard Poisson brackets 

{qm,qn} = i( l - < ^\qn\ 2 )§m,n, {^m, q n } = {q* m ,q* n } = 0. (4.8) 

The equations of motion are derived in the usual way, q m = {H, q m }. Using the follow- 
ing properties of Poisson brackets, 

{a,bc} =b{a,c} + {a,b}c (4.9) 
{f(a),b}=^{a,b}, (4.10) 



100 



the AL equation can be derived from the equations of motion of Hamiltonian (4.7), 



<?m —{H,q m } 

= (Unq* n+ l,qm} + {qtin+UVm}) ~ ^E{ ln f 1 ~~ ^nf) >*»} 

n ° n V Z 7 



=z 



M ® n 

(-0C 1 -?kn+l| 2 )8m,n+l + <?n+l • ("OC 1 - ^ l<?n| 2 )Sm,n 
n L z z 

4^ (o/2)q n a 

+ ^lT3^p-H)(i-2^l 2 ) 5 - 



{q* n ,q m } 



9m-l(l- ^kml )+<? m+ i(l--|<? m | ) 



~-i{qm-\ +qm+l — 2<? m ) — I^l^wl + <?m+l) 



This equation is completely integrable and has an infinite number of conserved 
quantities (for the infinite lattice) and can be solved by the method of Inverse Scatter- 
ing Transform (1ST) developed by Gardner, Greene, Kruskal and Mira (Gardner et al, 
1967). 



4.2 Conserved Quantities of the AL equation 

In this section, we outline the approach of inverse scattering and the work of Ablowitz 
and Ladik for calculating the conserved quantities of the AL equation (Ablowitz and 
Ladik, 1976; Ablowitz and Segur, 1981). 

The method of inverse scattering is analogous to the Discrete Fourier Transform. In 
the direct scattering problem, scattering data is derived from initial data, the potential. 
The time evolutions of the scattering data is simple. From the scattering data at some 
time 't', the potential can be calculated via the (non-trivial) inverse scattering transform. 
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In order to calculate the conserved quantities is it only necessary to calculate the scat- 
tering data. From the scattering data it is clear which quantities are time independent. 

The following is an outline of Ablowitz and Ladik (1976) with details filled in using 
(Ablowitz and Segur, 1981). For the AL equation, we taking S n (t) = T n (t) = in 
(Ablowitz and Ladik, 1976). The canonical pairs are R n and Q n . In the end we set 
R n = ±cc<2*. 

To begin, consider the generalized eigenvalue problem, 



Vi, n+ i =zV ljn + Q n (t)V 2 , n 
1 

V 2 ,„+i=-V 2 ,„ + R n {t)Vi, n , 

z 



(4.11) 



which can be equivalently expressed as 



E -Quit) 
-R n (t) E 





(4.12) 



where E is the shift operator: EX n = X n+ \. In this form, one can see the analogy with 
the Schrodinger equation in quantum mechanics with Q n and R n playing the role of the 
potential. The potentials correspond to the real-space classical fields \|/ w ,\|/£ in the AL. 
The time-dependence is postulated to have the form 



Vi, n =A n (t)V hn + B n (t)V 2 ,n 

(4.13) 

V 2 (t)Vi yn + D n (t)V 2yn . 
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The associated equations of motion for A„,5„,C„,D„ are generated by cross- 
differentiating (4.11) and (4.13) 

J t (Vl )B+ l) =ZV1,„ + Q n V 2 ,n + QnV 2 ,n 

=z(A n V hn + B n V 2 ,„) + QnV 2 ,n + Qn (C„V ljn + D n V 2 , n ) 
d , . 1 . 

57 {V 2 ,„+i) = -V 2 , n + R n V hn +R„V hn 
dt z 

= - (C n Vi jn + D n V 2 , n ) +R n V hn +R n (A n V hn +B n V 2 , n ) 

z 

dV ln >\ (4.14) 

=A„+iVi in +i+B n+ iV2,n+l 

=A n+1 (zVi, n + Q n V 2 ,n) +B n+1 ( -V 2; „ + R n Vi !n 



dt ' n'=n+l 



Z 



dV 2 „, 



dt 



= C n + 1 V\ tn+ 1 + D n+ 1 V 2;ll+ 1 

n'=n+l 



=C n+1 (zV hn + Q n V 2; „) +D n+l (^V 2/l +R n Vi^j 



where the eigenvalue, z is time-invariant and the explicit time-dependence has been 
dropped. Requiring f = (i^n') 



=n+l 



Z (A B V 1)B + fl M V 2 ,„) +QnV 2 , n + Qn (C B Vi )fI + £> B V 2 ,«) 



= A M+ i (zVi,„ + £„V 2 ,„) +B n+ i (-V 2!n + R n V hn 
! (4-15) 

- (C„Vi, n +D„V 2 ,«) +^nVl,n (A„Vi,„ +5„V 2) „) 

z 

= C„+l (zVl,„ + e„y 2 ,n) +D n+1 (^V 2 , n +R n Vl, n 
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and equating the coefficients of V^ n for each equation results in the following time- 
evolution equations, 



zA n A n +R n B n +i - Q n C n = 
(l/z)B n+i -zB n + Q n {A n+l -£>„) = Q n 

(4.16) 

zC„+i - ( 1 /z) C n - R n (A n - D n+ 1 ) = R n 
( 1 /z)A n D n + Q n C n+i - R n B n = 0, 

where A n X n = X n+ \ -X n . 

The linear dispersion relation and the above equations suggest the following expan- 
sions for the functions of (4.13): 

A n = z An ^ -\-A n ^ 

B n = zB n ^ + -B n "* 

(l) 1 (l) (4 ' 17) 
Cn = zC n H C n 

z 

D -Z)( ) + Iz)(- 2 ) 

L^n — LJ ii > n 

Z Z 
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The coefficients of these expansions are determined by subsituting back into (4.16), 

z 3 A n A {2) +zA n A i n ) +zR n B i n l l l + (l/z)R n B [ n l\ 

-zQnC { n ) -{\/z)QnC ( n i) =0 
^n+l + lVZ )B n+l -Z B n -B n +Q n {Z A n+l +A n+l 



-D^-(l/z 2 )D { - 2) )=Q n 
z 2 C% + C™ -<£> -(l/z 2 )d l) -R n (z 2 A^ +A^ 

-Dfl-{\/z 2 )D^)=R n 
(l/z)A n D^ + (l/z 3 )A,A 2) +zQ n d n l } 1 -(l/z)Q n d n l\ 

-zRnB i n ) -(l/z)R n B ( n l) =0, 



and solving in powers of V: 



a„a! 2) =o 




A (2) - 


-A (2) - 


a„d!" 2) = 


=> 


D (-2) 











tf > = ( 




= 


=> 


" — 


_n D (- 2 ) 






J*" = < 


r (-l) r, „(-2) . 


= 




> 4" = 



0(z): A^f = fi„d 1) -J?X+i 

= -AJ 2 HQ n+l R n -Q n R n -i) 
^A^ = -AS 2 )Q n R n ^+A_^ 
O ( : A n D^ = -/?„5i ^ — QrcC^j 

= -Dj- 2 )(e„7?„ +1 -e„_i J ??„) 



(4.18) 



(4.19) 
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The zeroeth order yields the time-dependence of the potentials, 

o (z") : 

G„= B < 1 » 1 - B <r 1) +&(Ai°» 1 -ri 0) ) 

Qn = (a_ (2) G„+i -D_ ( - 2 »e„-l) (1 -&■«„) + & (a_<°> -D_<°>) (4.20) 

R„ = (DWR„ + , -Aj- 2 >fl„_,) (1 - Q„R„) -R„ (a.I»> -D.l»>) . 
In summary, the coefficients of the time-evolution equations are given by 

A n =AS 2 Xz 2 -QnRn-l)+A_W 

B n =A_^Q nZ -DS^Q n ^(\/z) 

(4.21) 

C n =A_^R n ^z-DS^R n {\/z) 

D n = DS- 2 \(l/z 2 )-Q n -M+Dj°\ 

while the time evolutions of Q n and R n are governed by 

Qn = (A_WQ n+l -DS^Qn-^j (1 - Q„R„) + Qn (aJ°) (4.22a) 
R n = (pS 1 >R n+1 -Aj- 2 )/? n _i) (1 - (2„fl M ) - J? B (aJ°) (4.22b) 

4.2.1 AL equation 

To obtain the AL equation, we let 7?„ = ±a£>*, AJ 2 ) = -DS' 2 ^ = i, and A_(°) = 
— D_(°) = —i. In this case equation (4.22a) becomes 

2« = i [Qn+i + Qn-i- 2QnTa\Q n \ 2 (Q n+l + Q n -i)] , (4.23) 
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which is the AL equation. With these substitutions, the time-dependent function 
A n ,...,D n obey 



A n = i(z 2 ^aQ n Q* n _ l -l) 

B n = i{zQn-{l/z)Qn-l) 

(4.24) 

C n = ±ia(zQU-(l/z)Q* n ) 
D n = i{\-{\/z 2 )±aQ n ^Q* n ). 

4.2.2 Direct Scattering Problem 

Asymptotic Solution of Generalized Eigenvalue Problem First of all we assume 
that Q n and R n are on compact support, that is that as \n\ — > °°, Q n ,R n ^0. In the limit 
\n\ — > oo, the generalized eigenvalue problem becomes, 

{Vl } n+1 = zV\, n 
(4.25) 

To solve this direct scattering problem, define the time-independent eigenfunctions, 
(j)„,(j)„,\|/, 3 ,\j/„ which have the asymptotic forms: 



n — > —oa : 



n — > +oo : 




(4.26) 
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These eigenfunctions satisfy the generalized eigenvalue problem at a fixed time, which 
will be taken to be t=0. They do not solve the time-evolution equations, which we will 
return to later. One can establish by induction, that 

z~ n § m z"\|/ n are polynomial in powers of - and are analytic for Id > 1 

z 

z w (j) n , z~ n ^f n are polynomial in powers of z and are analytic for |z| < 1 (|z| < °°) 

(4.27) 

given that Q n and R n are on compact support. 

Wronskian and Linear Independence The Wronskian of a set of functions 
{<|>i, ...§ n } is defined by: 



w(4>i,...<M = 



4>i 



¥ 2 



4 n_1) 



(4.28) 



If the Wronskian is non-zero in some interval, then the functions are linearly indepen- 
dent in that interval. 
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For this problem, the Wronskian is defined as W n (u,v) = u\^ n V2,n ~~ W2,n v l,«- 

The 

Wronskian of the functions that obey (4. 1 1) can be found by 

W n+ i(u,v) = U\^ n+ \V2,n+\ -U2,n+\V\,n+l 

)((l/z) 

+ R n v l 

,n ) 

- {{\/z)u2, n +R n U\,n) {zv\, n + QnV2,n) 

= U^ n V2, n + zR n U\,nV\,n + (1 /z)Q n U2, n V2,n + Q n R n U2,nV\,n ( 4 '29) 

- ("2 + (\/z)QnU2,nV2,n + zR n Ul,nVl,n + QnR n U2,nV\,n) 
= {\-R n Qn){u\,nV2,n-U2, n V\,n) 

= (\-R n Qn)W n (u,v) 

Using the asymptotic forms of (|)„,^ n ,y M ,vj/„, one gets 

lim W n (M) = -h,n<$>l,n = -(-Z- n )(z n ) = 1 

(4.30) 

lim W„(\j/,\|/) = -\j/i,nV2,n = (z")(z" B ) = 1- 
By induction, one can show that the Wronskian of the eigenfunctions obey 

n oo 

W„(M)= -RiQi), W n (y,y)= J] (1-^ifii)" 1 (431) 

i=— co i=n+l 

for z on the unit circle. For i?„ = — Q*, W n is positive-definite so that <j>„, (|)„ are linearly 
independent. Otherwise assume that initially R n and Q n are less than one. Likewise one 
can show that vj/ M ,\|/„ are linearly independent. Thus one can define the scattering data 
by 

4>„ =a(z,t)y„ + b(z,t)y n 

(4.32) 

4>„ = -fl(z,?)Vn + ^(z,0¥n 
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where t is a parameter. 



Time Dependence The eigenfunctions satisfy the generalized eigenvalue 

problem, but not the time evolution equations (4.13). In the asymptotic limit (4.13) 
becomes 

Vl, n = A±Vl, n 
V 2 ,n = D ± V 2 ,n 



(4.33) 



where 



A± = lim (A n ) = z 2 Aj 2 ) +AJ°) 

D± = lim (D„) = (l/z 2 )D_ (2) + dS°\ 

n— >±co 



(4.34) 



Define a new set of eigenfunctions which will satisfy both the generalized eigenvalue 
problem and the time-dependence 



4>i f) = % exp(A J) \\f ( n ] = Mf n exp(D + 
ffi = <j>« exp(D_f) yi f) = \|/„ exp(A + ?) 



Both {tf } and {\|^ f ' ) ,\|/i f ' ) } are linearly independent, so we may write 



^ =a \j/P+M/? ) 

xW - (0 , r -W 



(4.35) 



(4.36) 



where ao,ao,bo,bo equal the scattering coefficients, a,a,b,b at t=0. Substituting (4.35) 
into (4.36) gives 

^=a e^+- A -^ + b eU D +- A -^ n 

(4.37) 

ft, = -doe^- D -^ n + b e^+- D -W% 



110 



Comparing with (4.32) and noting that A + = A_ and D + = D , the scattering coeffi- 
cients are 

a = a , b = b ett D +- A -W 

(4.38) 

a = a , b = boeK D +- A -W. 

The most significant result here, for our purposes, is that the coefficients a, a are con- 
stant. In the next section, we use this to derive the conserved quantities of the AL 
equation. 

4.2.3 Conservation Laws 

The conservation laws can be derived by considering the asymptotic form of a{z). 

% ~ -d(z) | ° j z~ n + b(z) (^\z n a(z) ~ -z n h,n (4.39) 

Next substitute § n into the eigenvalue problem (4. 1 1), 

z n ~%, n+ \ =z n h,n + z n - 1 Qnh,n (4.40a) 

Z n h,n+1 = Z n - l h,n+Z n Rnh,n. (4.40b) 

Solve (4.40b) for z w (j)i, n and substitute into (4.40a), to eliminate <j>i )W 

Z n h,n = l~ {Z n h,n+1 -Z n -%, n ) = -^A n (z>2,n) 
Kn ZK n 

f-2 J_ An ( Z »+^ 2;M+1 ) = -L A „ ( Z n h,n) +Z n - l Qnh,n (4.41) 
K n +\ ZK n 
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which can be written as, 



a ( n -2«-l/ «Z x 
A " ( z 2n+l^ ) - 2nZ (Z <t>2,n)- 



(4.42) 



Define 



-z n h, n = n 



(4.43) 



such that A M (-z>2,n) = A n (nL-oog*) = (gn+i - 1)I1L-ooS* and substitute into (4.42) 
to get the recursion relation 



gn+l(gn+2-l) -Z 2 %^(,?n+1 - 1) =Z 2 Rn+lQr, 



{AAA) 



which can be re-written as 



Rn-\gn{gn+\ - 1) -Z 2 R n {g n -\) = Z 2 R n R n 1 Q„- 1 . 



(4.45) 



Expand g M and in powers of z , 



g n = gn 0) +z 2 g ( n l) + z 4 g n 2) + ...= ^g n m) z 2m 

m=0 



(m) 2m 



gn+l = g [ n +l +Z 2 g„i l +Z 4 g [ „l l + •••=£ gT+\Z 

m=0 



(4.46) 



Substitute into the recursion relation for g n , g n +i, 



Rm £ gPz 2m (Z,g ( :Uz 21 - 1) -z 2 R n ( £ g n m) z 2m - 1) 

m=0 1=0 m=0 



= Z R n Rn \Qn-l 



Rn-1 



£ z 2(m+l) (m) a (l) 



m,l=0 



L Jim Am) 
Z gn 

m=0 



2 



7 2(m+l) W 



Z Z R n (Rn-lQn-l-l)- 



(4.47) 
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Solving recursively for g„ ) in orders of z 2 , gives 



Sn — 1 



0(Z 2 ): g { n ] =Rn-lQn-2 

0( Z 4 ) : ^ = J?„-lG„- 3 (1 -Rn-lQn-l) 

For orders of z 2p , p>2, this recursion relation reduces to 



Rn-l 



m=0 



A m )Ap- m ) _ Ap) 

in 6 n +l sn 



= R n - 1 



m=l 



or 




which gives expressions for higher orders of the coefficients 

0(z 6 ): 



p)_Rn-l o (2) _ (1) (2)_ (2) (1) 
>n-l 6„-lSn 8 n -i$n 



0(z*): 


* (4) 


_ Rn 


-1 


sn 


Rn 


-2 


0(z 10 ): 


£ (5) 
Sn 


_ Rn 


-1 


Rn 


-2 



*»-2 

Rn-\Qn-\{S - Rn-^Qn-3 - Rn-lQn-2 

+ Rn-zQn-^Rn-lQn-l) 
-Rn \Rn -iQl 3 (1 -Rn-lQn-l) 

n-1 <5„-l<5n S n -\sn S n -\Sn 



Next use the expansion 



00 a" 



log(l+x) =£-(-1) 



ra+1 



n=l 



to write out the series expansions, 



log[a(z)] = Umlog[ ]J g k ] = £ log(g„) = £ log 



£ gl m) z 2m 



m=0 



= £ log 



(m) 2m 



m= 1 



CO CO 



= E 1^ 

n=— <x>p—\ P 



p+l 



(m) 2m 



gn 

m=l 



p (4-53) 



Since a(z) is a constant of motion, each coefficients of z 2a in the above expansion of 
the must also be time-independent. These coefficients are the conserved quantities. 
Expanding the first few terms gives, 



CO CO 
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(4.54) 



^5(1 



+ o(z 12 ). 
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4.2.4 Conserved Quantities 



The first few conserved quantities are: 



c 2 = S>i 2) 4 



£ RnQn-2 ( 1 " Rn- 1 Qn- 1 ) + « ^fij- 1 



J 1 ) 



= (1 -Rn iQn 2 ~ R n -\Qn-\ + Rn-lQn-lRn-lQn-l) 
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— R n R n -\ Q Z n _ 2 (l -R n - X Q n - X ) +R Z n Q n -lQn-2(l ~ R n -lQn-l) 



+ \{RnQn-lf 



C4 = £^ 4) + ^ 



£ (2) 



<. (1) 



c 5 = Es< 5, +A» l +A, ( r , + 



e (2) 



^ (2) + - 
! £ (1 W 3) 



5« 



<5« 



+ 5 8n 



(4.55) 



Further conserved quantities can be determined by using the recursion relations and 
expansions given. From looking at these equations, we notice a few patterns. First of 
all, the leading order in a of C m is En and the largest order of a is a m . 

4.3 AL and BHM 

To obtain the AL, 4.23 we let R n = OLQ* n . Furthermore, to put it in a more familiar form, 
we write it in terms of the real-space fields, \j/ n = Q n , and vj/* = Q*. 



dt 



% = i [Of«+i +¥n-i -2\j/„) -a|\j/, 3 | 2 0q/„ + i + \j/„_i)] 



(4.56) 



115 



Recall that in real-space the BHM has the form 



—% = -i [~J (%+!+%-! - 2%) +IUQNs\%\ 2 ty] 



(4.57) 



By defining a new time scale, x = Jt, 



(\|/ n+1 +\|/ M _l -2\|/ n ) — |V«rV B 



(4.58) 



Note that here n is an index of real-space fields, following the notation of Ablowitz 
and Ladik (1976), while in earlier chapters it was used as a momentum index. From 
these two forms, it can be seen that these two equations are mathematically close for 
a = 2hqN s /J = 2kN s and identical in the limit \|/ n+ i ->■ 2\\f n . 

The first few coefficients of g n are 



Sn ] =a\j/*_ 1 \j/„_2 

gn ] =a\j/*_ 1 \j/„_ 3 (1 — cx|\j>„ 2 | 2 ) 

(3) (2 ) (1) (2) (2) (1) 

8n ~ ™V7V* °n-\ S n -\Sn S„-lSn 

= caj/*_ j \j/ M _ 4 ( 1 - a | 1 2 - a | \j/ M _ 2 1 2 + a 2 1 \j/„_ 3 1 2 1 vj/ n _ 2 1 2 ) 
-a 2 vj/^_i¥n-2¥L 3 (1 -oc|¥«-2| 2 ) • 



(4.59) 
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In terms of the real-space fields the first few conserved quantities are 



Ci =£a\fr*\fr„_i 



gn 



n Z 
<^3 — Lgn + gn gn + T 

11 



2 ^ 
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In terms of the momentum- space fields {\|/ p ,\|/*}, Ci and C 2 are 



C 1 = a£| Vp |V 27I! ^ 



= a£lYpl 



cos 



sin 
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+ 



2M 



(4.60) 



(4.61) 



As seen in earlier chapters, there is a threshold for chaos in the BHM. Below this 
threshold there are initial states that do not thermalize, but nevertheless relax to a steady- 
state. For fully integrable systems it is known that this steady-state can be described by 
a constrained thermodynamic ensemble that accounts for all of the integrals of motion. 
Anecdotal evidence suggests that the initial total quasi-momentum plays a role in the 
dynamics of the BHM for small nonlinearities, K. Is it possible that the other integrals of 



117 



motion of the AL affect the dynamics? Is the nonthermal steady-state governed by a con- 
stained ensemble that takes into account the conserved quantities of nearby integrable 
models? There are other nearby integrable models, which include the noninteracting 
case and the continuum limit. The imaginary part of the first conserved quantity of AL, 
Im C\ , is analogous to the total quasi-momentum and Re C\ to the total kinetic energy. In 
the noninteracting limit the momentum distribution is conserved, so any higher moments 
of the momentum distribution is also conserved. If the mapping between the BHM and 
AL fields corresponded to simply equating the fields (V|/„,bhm = V«,al)> then conser- 
vation of C\ in BHM would not distinguish between being close to the noninteracting 
case and close to AL. In contrast, the second and third terms of C2 are not conserved in 
the noninteracing model. However, the mapping between the BHM fields and the AL 
is nontrivial and the nonlinear corrections to the mapping are expected to allow one to 
distinguish between effects of quasi-conserved quantities of AL and of the noninteract- 
ing case. A derivation of the mapping between the BHM and AL fields and a study of 
the role of the conserved quantities of AL in the dynamics of BHM is the subject of 
proposed future work. 
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Chapter 5 



Outlook and Conlusion 

5.1 Summary of Results 

One of the fundamental assertions of statistical mechanics is that the time average of a 
physical observable is equivalent to the average over phase-space, with microcanonical 
measure. A system for which this is true is said to be ergodic and one can calculate 
dynamical properties of the system from static phase-space averages. Dynamics of a 
system which is fully integrable, that is has as many conserved quantities as degrees 
of freedom, is constrained to a reduced phase space and thus not ergodic, although it 
may relax to a modified equilibrium. What happens as one moves away from the fully 
integrable case? 

In this work we have studied the relationship between chaos and thermalization in 
the one-dimensional Bose-Hubbard model in the classical-field approximation. We have 
compared two quantitative measures of chaos and thermalization: (1) the finite-time 
maximal Lyapunov exponent, averaged over a microcanonical ensemble and (2) the 
normalized spectral entropy, which is a measure of equipartition of a modified energy 
in the independent mode approximation. There is a strong correspondence between the 
Lyapunov exponents and normalized spectral entropy. 

We find a threshold for chaos and a corresponding broad transition from incomplete 
to complete thermalization. The stochasticity threshold is governed by two parameters: 
the strength of the nonlinearity K and the average energy-per-particle, e T . Both of these 
parameters are finite in the thermodynamic limit and suggest that the threshold will 
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survive in that limit. Far above the threshold, in the strongly chaotic regime, relaxation 
to the thermal state is complete. In this region, the fluctuations in kinetic energy scale 
as y/N~ s confirming their thermal nature. We study the size scaling of the Lyapunov 
exponent and find that it is universal with respect to the size of the lattice. For small 
nonlinearities, the stochasticity and thermalization thresholds are finite for the range of 
energies studied and don't show tendencies to vanish at high energies. 

In the vicinity of the threshold the relationship between chaos and thermalization 
is complex. There is a transient regime supporing both chaotic and regular trajecto- 
ries, so that the Lyapunov exponent is non-zero, but full thermalization does not occur. 
Remarkably, in this region individual initial states with larger Lyapunov exponent tend 
to relax closer to the thermal state. There is also a region where although the Lyapunov 
exponent is zero, there is significant relaxation towards the thermal state. We conjecture 
that this redistribution in the phase-space is due to the quasi-regular dynamics governed 
by the nearby Ablowitz-Ladik lattice, which is integrable. Above £7- ~ 0.67, both the 
stochasticity threshold and thermalization threshold overlap very closely and appear to 
depend only on the nonlinearity strength. This suggests that there is a critical nonlinear- 
ity, K c . ~ 0.2, such that the system is regular for K < K c independent of the energy of the 
system. In the opposite limit of small energy and large nonlinearity, there is a separation 
of thresholds, where almost complete relaxation to the thermal distribution is observed 
in the absence of chaos. 

An analysis of resonances of the BHM gives a Chirikov's criterion for the chaos 
threshold that depends on parameters which vanish in the thermodynamic limit. This 
conclusion is confirmed on dimensional grounds. The criterion predicted by the Chrikov 
criterion is different from the one inferred from numerical calculations, signifying the 
failure of the standard Chirikov's approach. 
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Anecdotal evidence suggests that the total quasi-momentum may play a role in the 
relaxation dynamics. The quasi-momentum is strictly conserved in the nearby fully inte- 
grable Ablowitz-Ladik model, as well as in the non-interacting and continuum limits. 

There are at least three known near-by integrable models: the Ablowitz-Ladik lat- 
tice, the continuous nonlinear Schrodinger equation and the noninteracting model. We 
outline the method of Inverse Scattering Transform and generate all of the integrals of 
motion of the closely related, fully integrable model of Ablowitz-Ladik. We conjec- 
ture that the presense of quasi-conserved quantities may alter the scaling of the chaos 
criterion. 

5.2 Open Questions 

These observations lead to many questions that deserve further investigation. 

• What is the reason for the failure of the Chirikov criterion to accurately predict the 
chaos threshold? Is it related to interference between resonances due to near-by 
integrable systems? Could the proper scaling be recovered in a multiple-resonance 
model? 

• What is the underlying theory that governs the threshold? 

• For K > 1, how does the chaos threshold scale in the thermodynamic limit? Is the 
number of modes involved relevant? 

• What governs the slow relaxation times where they appear? Is it related to the 
conserved quantities of the near-by integrable systems? If this is the case, which 
near-by integrable system? 
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• For those states that show no signs of thermalization, what governs the steady- 
state? Can these states be described by a constrained ensemble? Which "quasi- 
conserved quantities" are the relevant for the constrained ensemble? 
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Appendix A 



Thermodynamic Distribution within 
Hartree-Fock 

A.l Hartree Fock 

In this section, the following integrals are used: 

Jo 



dxe-™ = - 
a 

1 



L 



d xxe - ax = ^ (A.l) 
o a 2 

dxx 2 e- ax = 4 
o a 3 



Note that N is the number of degrees of freedom and a = L/N is the lattice spacing. 

Within Hartree-Fock the form of the density distribution function is taken to be 
Gaussian and the thermal expectation value of the Grand Potential, 

(F) = {H)-T{S)-v{N a ), (A.2) 

is minimized, where N a is the norm. The density distribution function with two-body 
interactions has the form, 

a H F = ^ exp ^£ -a n y\\f n \\f* n }j = ^ exp (- L a « 7 «) • ( A - 3 ) 
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In the second step, we assume that off-diagonal terms are zero. The a n coefficients are 
unknown and are determined by the condition of minimizing the grand potential. The 
partition function Z, normalizes Ohf so that the integration of Ohf over all of phase 
space is 1. Throughout the sums run from 1 to N, where N is the number of lattice sites. 



Z=-J— f d N B f d N Ie-^ a "'" 
(2nh) N J J 

j r27Z p2tZ r27Z r°° r°° poo 



1 poo poo poo 

(A.4) 

=wL dhe ~ aih L di * aih -L d! ^" 

N l 

The expectation values of each term in the Grand Potential is calculated, using the 
Hartree-Fock density distribution function, Ohf- The expectation value of a generic 
observable is given by 

{0) Z d% l dNl0 ^n})OHF 

= J^t\ ai l dNe l dNl0e ' LnanIn 

The expectation values of the relevant observable are calculated below. 



(A.5) 
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Norm 
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Hamiltonian - Kinetic Term 



i n r r \ 
(Ho) =J^EfU ha iJ dNQ J dNl L 1 ^ 

N (£> 

w n 

n=\ a " 



Hamiltonian - Interaction Term 
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Entropy 

^--^I** I ********* 
=^flm I * N e- M (-EOA + Iogn^) (A.7) 

n a « j 

=N-^log(haj) 

j 

A.2 Minimization of the Grand Potential 

The thermal expectation value of the Grand Potential within Hartree-Fock is given by 

m II m p vx m KX p \ m / m UUjm 

Taking the variation with respect to a n , and setting it equal to zero gives 

m = .^L 2L ± + L + jL =0 (A . 9) 

8a„ a 2 . h 2 al^a m a„ ha 2 
Using £ m a m ' = hN a and solving for a n , 

a " = -j^ l h(S >n + 2 VoN a - n\ (A. 10) 

The thermal expectation values of the occupation of momentum mode n become 

1 Pi T 

//„) = _ = . (A.11) 
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In general, the coefficients fi and T are unknown and are determined by imposing con- 
straints on the norm and energy, which come from the dynamical code. The constraints 
are 

n 

(A. 12) 

Et = (H) = £ co„ (I n ) + g £ (7 m > (/„> = £ to„ (/„> + m N 2 a 
Beginning with the expression for (/„), we can solve for T in terms of jji, N a and energy. 



T = co M (/„) + 2^N a (I n ) -^(I n ) (A. 13) 



Summing over n 



T=^[E k + 2 l j N 2 a -nN a ] {AAA) 



where E k = (H n (I n ). This expression for T can be substituted back into the constraints 
to reduce the system to two equations with two unknowns. Using the expression for 
temperature and normalization condition, a single constraint remains to be solved, 

1 [E k + 2wN*-yN a ] 

[fm n +2^N a - M ] Na ~° (A - 15) 
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